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We consider a single-component gas of dipolar bosons confined in a one-dimensional optical lattice, 
where the dipoles are aligned such that the long-ranged dipolar interactions are maximally repulsive. 
In the limit of zero inter-site hopping and sufHciently large on-site interaction, the phase diagram is a 
complete devil's staircase for filling fractions between and 1: every commensurate state at a rational 
filling is stable over a finite interval in chemical potential, and for every chemical potential the 
system is in a gapped commensurate phase. We perturb away from this limit in two experimentally 
motivated directions involving the addition of hopping and a reduction of the onsite interaction. The 
addition of hopping alone yields a phase diagram, which we compute in perturbation theory in the 
hopping, where the commensurate Mott phases now compete with the superfluid. We capture the 
physics of the Mott-superfluid phase transitions via bosonization. In a finite trap, we argue using 
an LDA and simulated annealing that this results in regions of commensurate states separated by 
patches with dominant superfluid correlations. Further softening of the onsite interaction yields 
alternative commensurate states with double occupancies which can form a devil's staircase of their 
own; we describe these states and discuss the possible transitions between them. Adding a hopping 
term in this case produces one-dimensional "supersolids" which simultaneously exhibit discrete 
broken symmetries and superfluidity. The contents of this work was first published as a chapter 
in the doctoral thesis 'On Exotic Orders in Strongly Correlated Systems' (Princeton University, 
09/09, supervised by S. L. Sondhi), and constitutes a considerably more detailed version of [PRB 
80: 174519 (2009)]. 
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I. INTRODUCTION 



When commensuration effects compete with long-ranged interactions, a starthng richness of phases can arise. This 
interplay between long-ranged forces and lattice effects has been studied in a variety of models - from the Ising model 
[H to interfaces between crystal surfaces Q ■ In this work, we explore a new context in which such physics leads to rich 
mathematical structure ~ that of the cold dipolar Bose gas. Recent advances in laser cooling and trapping have opened 
the possibility of creating cold Bose gases in optical lattice potentials. If these bosons have sufSciently strong dipolar 
interactions, the phase portrait is controlled by the interplay between the infinite-ranged dipolar repulsionQ, and 
commensuration effects due to the optical lattice - giving a potential experimental realization of the diverse phases 
expected in such a system. Furthermore, considering which parameters can be easily tuned in such experiments 
naturally opens questions about the ground states of these systems in new regions of parameter space. 

The driving force behind the strikingly rich phase diagram in these systems is that long ranged interactions on 
a lattice can stabilize exceptionally intricate ground-state structures in the classical, or strong coupling, limit. A 
particularly elegant example of this is the case of classical particles in a one-dimensional (ID] lattice interacting via 
an infinite-ranged convex potential studied by Pokrovsky and Uimin, and Hubbard (PUH) |3,@1- Here, it can be 
shown that the ground state filling fraction as a function of chemical potential /j, is a complete devil's staircase @, in 
which every rational filling fraction between and 1 is stable over a finite interval in fi, and the total measure of all 
such intervals exhausts the full range of fi. 

This devil's staircase of PUH has dramatic consequences for the physics of quasi ID cold atomic gases. Building 
on the existing understanding of this classical limit, we consider two perturbations of the devil's staircase that arise 
naturally in the experimental setting of cold atomic gases. The first of these is the introduction of a quantum kinetic 
energy, due to the finite depth of the optical lattice, which now renders the problem sensitive to particle statistics. 
We focus on the bosonic case, as this case is most likely to be realized experimentally. Tuning away from the classical 
limit has the well understood effect of initiating a competition between the crystalline, Mott phase that exists at zero 
hopping and the superfluid (Luttinger liquid in d = I) that must exist at all fillings at sufficiently large hopping. 
The phase diagram we find is thus an extremely complex variant of the usual Hubbard model story, with an infinite 
number of Mott-Hubbard lobes corresponding to commensurate phases. 

The second perturbation involves tuning the onsite interaction independently of the dipolar potential. This allows 
for a controlled departure from convexity, and generates new classical states comparable in complexity to the PUH 
states considered previously. This introduces doubly occupied sites in the classical limit. While describing the resulting 
phase diagram in complete and rigorous detail is beyond the scope of this project, we give an account of the "staircase" 
structure of the initial instability and of regions of the phase diagram where the classical states exhibit superlatticcs 
of added charge built on underlying PUH states. At least some of these regions exhibit devil's staircases of their own. 
Finally, upon introducing hopping we are led to an infinite set of "supersolids" — which in this context are phases that 
are both Luttinger liquids and break discrete translational symmetries. 

The structure of this paper is as follows. We begin in Section |lT] with an overview of the model and a sketch of 
the relevant experimental parameters. TH] gives a detailed description of the intriguing infinite family of classical 
solutions to this model, described in [l|, Q, and discuss the fractal structure known as the devil's staircase which 
constitutes the classical phase portrait. Then in Section [IVI we use a strong coupling perturbation theory similar to 
that previously studied in the Bose- Hubbard model , and in extended Bose- Hubbard models with nearest- neighbor 
interactions [1, , to calculate the boundaries of the Mott-Hubbard lobes where these commensurate classical states 
liquify. We also review how bosonization techniques admit a Luttinger liquid description of these phase transitions. 
This approach proves helpful in understanding the expected density profiles of atoms in a parabolic trapping potential 
- a question we address for both quantum and classical models in the convex limit in Section [V] 

In Section [VIl we switch gears and evaluate the classical phase diagram as Uq is tuned away from convexity. This 
reveals a mind-bogglingly complex array of possible phases, only some of which we are able to describe here. However, 
one class of such phases, in which the ground states resemble a two-component version of the HUP ground states, 
are of particular interest. Section IVlIl discusses in more detail the classical structure of these phases, as well as their 
behavior at finite hopping - where we find super-solid like phases. We have relegated several detailed calculations to 
Section [IX[ which contains supplementary material. 



II. THE ULTRA-COLD DIPOLAR BOSE GAS 



The unprecedented control over experimental parameters in trapped ultracold atomic gases has substantially 
widened the range of possible exotic phases of matter that can be explored. Optical lattices can be used to sim- 
ulate simple lattice models and/or vary the system dimensionality, while the interatomic interactions can be varied 
via a magnetically-tunable Feshbach resonance [lo| . Thus far, the focus in these systems has largely been on contact 
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interactions, since these generally provide a good description of atom-atom scattering in the low energy limit. How- 



ever, with the realization of atomic gases with strong magnetic dipole moments jllj and the prospect of working with 
molecules that exhibit electric dipole moments [H, , there is now substantial interest in examining the physics 
of long-ranged interactions in these systems [l^l- In this section we will present the model, motivated by these new 
experimental possibilities, which we will discuss in the remaining sections of this paper. 



A. Bosons in ID optical lattices 



We begin by discussing the basic physics of dipolar molecules in an optical lattice. Experimentally, the scenario is 
as follows: atoms, or in some cases even molecules, can be trapped using spatially varying electric fields. In optical 
traps, this frequency can be tuned to be close to an atomic transition. In this case the atoms are highly polarizablc at 
this frequency, and hence feel a relatively strong trapping potential. A lattice potential can then be made for the gas 
of atoms, by using counter-propagating laser beams to create a periodic pattern of maxima and minima of the electric 
field's intensity. The depth of the lattice is tuned by the intensity of the laser light; the lattice spacing is determined 
by the wavelength of the standing light-waves. 

To explore ID physics, in practice a 3D optical lattice is generated, but with a lattice depth and spacing much 
greater in the first two dimensions than in the third. This effectively creates an array of ID tubes, with a relatively 
weaker periodic potential and shorter inter-particle spacing along each tube. Typically, the whole ensemble of tubes 
is also confined in a long parabolic trap, to ensure that the particles remain inside the optical lattice system. 

Here we will ignore the complications due to interactions between tubes. To simplify the discussion, we will also 
begin by omitting the parabolic trapping potential - we will return to its effects in Section |Vl Hence the Hamiltonian 
we study is: 

i<j '^3 i 

-fj.'^rii - t'^clci+i + h.c. (1) 

i i 



This model describes bosons in a deep ID optical lattice with hopping amplitude t, on-site interaction energy U, and 
an infinite-ranged dipole-dipole interaction Vo/r^. Since the hopping potential t is controlled by the depth of the 
optical lattice, it can be tuned over a wide range of values. In the cases where the bosons are atoms, e.g. ^^Cr [ll[, 
the on-site interaction U may also be easily tuned using Feshbach resonances. We set the dipolar interaction to 
be maximally repulsive by aligning the dipoles perpendicular to the ID chain. We also note that though dipolar 
interactions will couple the different tubes, aligning the dipoles at an angle such that cos{9) = l/VS cancels this 
interaction in one of the two remaining dimensions. Hence to access the ID regime, one would ideally work with a 
single 2D array of optical tubes. 

Though our focus here is on the realm of possibly, rather than currently, attainable states of the dipolar Bose gas, 
it is instructive to ask what range of parameters can be attained in current experiments. The critical ratios are t/Vo 
and Uo/Vq. 

In cold atomic gases, there are essentially two ways to generate dipolar interactions: by cooling atoms, such as 
Cr, with large magnetic dipole moments, or by creating molecules with electric dipole moments - typically by 
generating heteronuclear molecules which have an intrinsic dipole moment, for example, with ^^K ^^Rb. Techniques 
for cooling and trapping atoms with large magnetic dipole moments have already been developed [isj ; cooling dipolar 
molecules is significantly more difficult, though significant progress towards creating cold molecules in their 2-body 
ground states has been achieved, most notably in fermionic systems [l6l . [l7| . Since the experimentally attainable 
electric dipole moments are much larger than their magnetic counterparts, it is interesting to consider the parameter 
r'egimes accessible in both types of experiments. 

Current ultracold bosonic molecules have electric dipole moments of the order of d « 1 Debye; for ^^Cr and ®''i?6 
the magnetic dipole moments are 6/iB and l/i_B, respectively. Using these values, the dipolar interaction strength at 
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a distance of n lattice spacings is 

_ 1 _ ^2^3.338 X m-^"Cmf 



1 



1.0 X IQ-^^Nm 



V'KIaq 



,,2 

« -^8.6 X lO-^^A^m (2) 

where a is the optical lattice spacing in nm, d = 0.6 is the electric dipole moment measured in Debye, and /x = 6.0 
for ^^Cr and 1.0 for ^'^ Rh is the magnetic dipole moment in Bohr magnetons. 

To find the hopping parameter requires calculating the overlap integrals of Wannier functions as a function of the 
depth of the optical lattice. The result is that the hopping t is given by [isj : 



4i?. / A ^ 



e 



where A is the intensity of the laser beam, and Er is the recoil energy of a single atom, Er = h'^/{8ma^). Current 
experiments can attain lattice spacings of order a = 500nm, at laser intensities on the order of 10 — 20 recoils. This 
gives a range of hopping parameters listed in Table \n~K[ We see that, assuming that the electric dipole moments of 
the molecules do not significantly alter the facility with which they can be trapped (a reasonable assumption since 
the frequency of the molecular transitions is significantly lower than that of the atomic transitions to which the laser 
light is tuned (l9|). fractions t/Vb 10"^ are well within the range of experiments on polar molecules. For magnetic 
dipole moments, the numerical value of t/Vo is larger by a factor of 300 (for ^^Cr) to 6000 (for ^''Rb), and hence we 
do not expect to see any Mott physics associated with these states in currently realistic trapping potentials. 
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300 


400 


500 


A/Er 








5 


0.184 


0.25 


0.306 


10 


0.048 


0.065 


0.081 


15 


0.016 


0.021 


0.026 


20 


0.006 


0.008 


0.010 



TABLE I. Experimentally attainable values of the hopping t/Vo as a function of the lattice constant and laser intensity in the 
approximate range of experimentally realizable values. Here we use parameters relevant to polar molecules - an electric dipole 
moment of 0.6 Debye, and the mass of '''^K *^Rb. The corresponding values for ^^Cr are larger by a factor of « 300. 



The effect of deforming the on-site repulsion U away from convexity is rather more difficult to probe experimentally. 
To attain this regime, Feshbach resonances must be used to tune the scattering length to be relatively large and 
negative, to compensate for the strong dipole-dipole interaction felt by two bosons confined to the same lattice site. 
Calculating the overlap of the Wannier wave functions in a deep optical lattice, where the particles are effectively 
localized at a single lattice site, gives 

U=^eJ^Y' «5xlO-2«^iVm (4) 
a \Er J a 

where is the scattering length of the bosons, tuned by Feshbach resonances, and we take = 20, a = 500nm. To 
this we must add the effective dipolar interaction, whose approximate order of magnitude we obtain by noting that 
two bosons confined to the same site are a separation of no more than a/2 apart - giving a dipolar energy of order 
2 X 10~^'^Nm. Hence a priori one would expect that on-site interactions are easily tuned away from convexity by 
making the scattering length sufficiently negative. In an optical lattice a reasonable limit [lO[ is as < a, so some 
tuning away from convexity may be possible[20j. 
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III. CLASSICAL SOLUTIONS AND THE DEVIL'S STAIRCASE 

The notion that longer-ranged interactions can stabihze fractional fillings in the Mott state is well-known in the 
study of extended Hubbard models. With only on-site interactions, the Hubbard model itself has commensurate states 
only at integer fillings; extended Hubbard models with repulsion between neighboring sites can stabilize such states at 
half filling, and so on. The HUP case is simply the extreme limit of the extended Hubbard model: convex interactions 
favor arrangements in which particles are spread out as homogeneously as possible given the filling fraction. In 
the case of infinite-ranged convex interactions, this yields a constraint that must be satisfied at all distance scales, 
leading to a pattern that is unique, up to global lattice transformations, for every filling fraction. Hence the classical 
configurations are devoid of all local degeneracies, and commensurate structure exists at all length scales. 

In our model Hamiltonian ([T]) , if [/ is sufficiently large the potential is everywhere convex - meaning that 

V{x) < XV{x - I + A) + (1 - X)V{x + A) for < A < I . 

The t = (classical) ground states of ([T]) in this case are therefore those of a PUH Hamiltonian: for every rational 
filling fraction ly = p/q, the ground state is periodic with period q [l[ . Each such ground state is unique up to global 
lattice translations j2ll |. We denote these states commensurate ground states (CGS). Adding or removing a single 
particle from a CGS in the infinite volume limit produces a q-soliton state (qSS) containing q fractionally charged 
solitons of charge 1/q. We now review the nature and energetics of the CGS and qSS states. 

A. Commensurate Ground States 

We now review the detailed construction of the CGS classical ground states. To do so, it is convenient to characterize 
a state by the set of all of its inter-particle distances. Call particles /*'' neighbors if there are I — 1 occupied sites 
lying between them (their separation will generally exceed this). Any static configuration of particles on the lattice 
is characterized by the set {S) of sets 5"; = {rp^rp-*, ...r\^^} of l*^ neighbor distances between a particle and its l*^ 
neighbor to the right. (On a finite-sized chain it is simplest to consider periodic boundary conditions, in which case 
rf' exists for all i and I < N.) The energy of such a configuration is 

I i 

where V{r) is the (convex) potential. Here the second sum is a sum over individual particles; the first sum is a sum 
over all neighbors of each particle. A solution which minimizes the inner sum for each I also minimizes the energy; 
hence it suffices to find the set Si for which 

Ei^J2^irl'^) (6) 

i 

is minimal, and establish that a configuration for which Si has this form for all I does indeed exist. 

1. Energetics of convex interactions 

The first step in constructing the HUP ground states is to establish which Si, for a given Z, minimizes the energy 
(pl). For convex interactions, the optimal Si is a maximally compact set - e.g. the one for which all r/ arc as close as 
possible to the same value. This is commensurate with the intuition that convex interactions favor maximally uniform 
distributions of particles, and we will prove presently why this is so. 

For a given I and filling fraction v, what is the maximally compact 5/? In the absence of a lattice, maximal 
compactness would imply r^*' = n for all i; however on a lattice we must also allow Si = {ri,ri + 1}. To see why, 
recall that the sum of all distances between first neighbors must be the total length of the lattice - or, more generally, 
that for a lattice of length L we have 

(7) 

i 

The average separation between l^^ neighbors is thus: 

^"' = ^E-i^^=V- (8) 
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where V = N / L. For rational fillings v = p/q, this gives: 

fl = Iq/p (9) 

Thus either p\l and the maximally compaet set is Si = {fi\ = ^g, or r/ is not an integer, and hence Si must have 

at least two elements - r'f^ € {r;,r/ + 1}. Hubbard [l| showed by construction that for all rational there exists a 
solution in which Si = {r/, r; + 1}, and hence that these are ground states. 

Before outlining the details of this construction, let us show why maximally compact sets Si correspond to minimal 
energy solutions. For convex potentials V{r) 

V{ri +x + y) + V{ri - x) > V{ri) + V{ri + y) (10) 

In other words, given 2 pairs of points with the same average, the average of any convex function V will be smaller 
over the pair which is closest together. This is simply a reiteration of the familiar definition of convexity for continuous 
functions, stating that the line joining any 2 points on the graph of V{x) lies above the graph. 

Now, given two distinct Si with the same average n, the convexity of the potential (jlOp ensures that the set with 
the narrowest distribution of its elements will have the minimum energy. More precisely, let Si be the maximally 
compact set at level I, containing either one or 2 elements, and let 

^ n = Mn + Nin + 1) . (11) 

Any other set Si with the same mean must by definition contain at least one element rp"* such that either r|^'' > r; + 1 

or r[^'' < ri. (If there is only one element in Si then of course both limits are given by ri). Let us assume that the 

former case holds, and that there exists an 7'^' = ri + 1 + x for some integer x. As both sets have the same average, 
this means that 



^ ri=Mri + {N -l){n + l)~x (12) 

In other words, the mean of the remaining elements must be shifted correspondingly to the left. Picking one of the 
shifted elements rp'' < r/ (where strict inequality holds if S* contains only r;), create S^^'^ by shifting rp-* -> r-j^-* — 1, 

(2) (2) 

— > + 1 and leaving all other elements of Si unchanged. Then 

EiSi) Eisl'^) = Virl'^) + F(rP)) - Virl'^ - 1) - V^(rf ' + 1 > (13) 

where the inequality follows from the convexity of the potential, since for rp'' > rp'' the elements of sj^^ are more 
narrowly distributed about the mean than the elements of Si. This construction can be repeated so long as some 
rp-* > ri + 1 exists - in other words, so long as Si ^ at each step in the construction the new configuration has 
lower energy. Since this holds for all I, the optimal solution is characterized by S* = {5*;}; with each Si maximally 
compact. 



2. Constructing the CGS 

Hubbard [l| gives a construction which, for any rational filling p/q, constructs a configuration of particles on the 
lattice for which every set of neighbor distances is maximally compact. A thorough discussion of this solution is 
given in [2l| ; here we present an overview. 

Hubbard's solutions are periodic, and can be expressed in the form: 

nin2...nj (14) 

where rii represents the spacing between the i"* and {i + iy* occupied sites. Given a continued fraction {no, rii, Uk}, 
let 



^ = l/(no + l/(ni + ... + l/n,)) 



(15) 
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Filling 


CGS 


particle soliton 


hole soliton 


1/9 


q 


q-1 


9 + 1 


2/5 


23 


22 


33 


3/7 


223 


222 


323 


5/12 


23223 


22322 


32323 



TABLE II. Occupancy patterns in a few CGS 



be the numerator and denominator of the related continued fraction {no, rii, ■■■rii}^ i < k. Now consider 

Xa = HQ 

Yo=no + l (16) 

with 

Xi = (Xi_i)"' ^Yi^i 

Y, = (X,_i)"'K,_i (17) 

for all i > 1. The Hubbard solution is given by Xk- A few examples of these CGS states and their solitons are given 
in Table IIIIA2] 

One can show by induction that this solution has the correct filling fraction. The numerator and denominator of 
the full continued fraction can be expressed in terms of numerators and denominators of shorter continued fractions, 
via [H 

Pk = nkPk-i +Pk-2 

qk ^ nkQk-i + qk-2 (18) 

Using ()16|) . we have Xq ~ no which contains one particle per no sites; Yq ~ hq + 1 contains one particle per no + 1 sites. 
Equation (jl7p states that if Xi-i contains Pi-i particles in g^-i lattice sites, and li-i contains Pi-i + Pi-2 particles 
in + qi-2 lattice sites (with p_i = 0, q-i = 1), then Xi contains riiPi-i +Pi-2 ~ Pi particles in Uiqi-i + qi-2 = qi 
lattice sites. Yi contains (n; + l)pi-i +Pi-2 = Pi +Pi-i particles in {m + l)qi-i + qi-2 = qi + qi-i sites. Thus at every 
step in the process the solution Xi has filling fraction ^ corresponding to the truncation of the continued fraction 
expansion at n^. 

Further, it is easy to convince oneself by inspection that these solutions do indeed satisfy Hubbard's criterion that 



B. The Devil's staircase 



We have seen how Hubbard's construction gives the ground states of any infinite-ranged convex potential at rational 
filling fractions. We now address the issue of these states' stability. Bak and Bruinsma showed that the range of 
fi over which each CGS is stable is given by 

V + - ^ (19) 

{nq + 1)3 (nq - 1)3 (nqf ' ^ ' 

At a filling fraction of ^ the result is independent of p, and falls off sharply as a function of g. At t = these intervals 
cover the entire range of /i pertinent to fillings less than unity, giving the devil's staircase structure. 

To understand Eq. (fT9l one must consider the effect of adding or removing a single particle from the CGS. We shall 
see that at filling fraction ^ this produces q fractionally charged solitons of charge l/q, each of which distorts the 
periodic ground state by altering the distance between one pair of adjacent particles by 1. (For every commensurate 
state, there is a unique distortion of this type which minimizes the potential energy). This results in a q-soliton state 
(qSS). 
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C. Structure of the q-Soliton State 

To understand why charge fractionahzes, consider the possible Sp after one charge is added or removed from a chain 
of length L containing particles. At rational filling N/L = p/q, for I ^ mod p. we have 

n., e{r,,n + l} (20) 

for all Z*'' neighbors i and j in the HUP CGS. For any value of /, the CGS must satisfy 

f^r, + A^n+i = N 

Nr,ri + Nr,+iiri + I) = IL (21) 

where TV^, , -/Vrj+i are the number of l^'^ neighbor pairs separated by distance r; and n + 1, respectively. After adding 
or removing a particle, the most energetically favorable lattice configuration requires Si be maximally compact for 
every I. For I < p, this requires that the set of possible radii {ri,ri + 1} remain unchanged; for I = p, however, we 
must now assume Sp = {q,q±l}. These newp*'* neighbor separations constitute the solitons of the qSS. After adding 
a particle, 

K + K+i = ^ + 1 

Kri+N^^^,{ri + l) = lL (22) 
giving the net change in the distribution of Z*'' neighbor distances: 

K - = ±(n + 1) 

Nr,^. - K+i = ±n (23) 

where +,— correspond to adding or removing a particle, respectively. Substituting rp ~ q (one charge removed) or 
Tp = q — 1 (one charge added) into ([231) shows that the qSS contains exactly q solitons. Thus each soliton has a 
positive (or negative, in the case of holes) charge of l/g with respect to the parent lattice. 

In a finite system, the qSS is also a HUP state with rational filling, obtained by distributing the q solitons as evenly 
as possible on the chain in order to satisfy Hubbard's criterion at all I. The denominator of such a state will be of the 
order of the number of lattice sites. It is useful to conceptualize filling fractions whose denominators are comparable 
to the system size as being qSS states of a related CGS of smaller denominator. 

Energetically speaking, in a finite system one must also account for the repulsion between solitons; for most of what 
follows, we will drop these terms and consider the infinite volume limit in which solitons arc infinitely far apart. In 
this limit the HUP qSS consists of q free solitons in the 'background' lattice of the parent HUP state. 



D. Energetics of the qSS 

We now turn to the question of the energetics: over what range of chemical potential is each CGS solution stable 
against the formation of solitons? Eq. above shows that the energy costs of adding and removing particles are, 
respectively, 

E+ = -fi+ [{ri + l)V{n)-riV{ri + l)] 

;7to(mod p) 

+ ^ nqV{nq - 1) - {nq - l)V{nq) + ... 

l—nq 

+ ^ ngy(ng+ 1) - (ng + l)y(ng) + ... (24) 

l—nq 

where we ignore soliton repulsion terms present in a finite system. Setting E± = in Eq. I24[ we obtain the values 
fJ.L{p/q) and fJ.R{p/q) of the chemical potential at the left and right extremities of the plateau. The width of the 
plateau (Eq. I19p is given by the difference between these. The sum over I ^ 0{modp) contributes equal amounts to /i^ 
and fiR, since Z*'* neighbors separated by r; and r; + 1 both exist in the initial solution, and creating solitons simply 
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adjusts their relative frequencies. However, the sohtons replace some Tp by rp ± 1 - giving different contributions to 
particle and hole like solitons. This gives the width of the plateau quoted in Eq. (|19p : 

MpIi) - t^Lip/q) = E [^("■'^ - 1) + ^("9 + 1) - (25) 

n 

which, by the assumption of convexity, is strictly positive. Since r„p — nq, the range of stability depends only on q. 

On the infinite-length chain, the function v^i) has a complete Devil's staircase structure. This means firstly that 
the function is monotonic and contains no finite jump discontinuities [2l|. Secondly, the set of all such intervals for 
rational fillings p/q with q > I completely covers the interval < /i < E-{v = 1), and in the infinite chain limit all 
ground states are periodic HUP states at rational fillings. 



E. Proof of devil's staircase structure 



Here we examine a few of the details required to show the devil's staircase structure. The fact that i>{ii) has no jump 
discontinuities essentially follows from the fact that between any pair of rationals on the real line, there is another 
rational. We thus focus on the second claim - that all ground states between v = Q and v = 1 are commensurate 
(except possibly on a set of measure 0). To do so, it suffices to show that the intervals of stability are disjoint, and 
that the sum of their lengths is the length of the relevant interval in ^. 

To show that the ranges of ji over which each rational filling is stable must be disjoint, we first point out that as the 
potential ^ is a convex function of the inter-particle spacing, the energy is a convex function of the filling fraction. 
Indeed, were this not so, it would be energetically favorable at some filling fraction to break a system of length L at 
filling ^ into(say) two subsystems of filling ^ + a; and ^ — x ~ which, as we know from Hubbard's solution, does not 
occur. This is because if e(p) is the energy per unit length of the state at filling p, then: 

2e(^) < e(^+x) + e(^ -x) (26) 

and the energy is a convex function of filling fraction. 

For any p'/g' < p/q^ choose I such that p' /q' = p/q — {l + l)/L, so that the state at filling [N — l)/L contains ^ji-l-l 
particles. The upper boundary of stability of the state p' /q' then occurs at the chemical potential where the energy 
of this state is 0. Let E{N/L) denote the lattice contribution to the energy of the state at filling N/L, such that the 
total energy of this state is iJ,N + E{N/L). We have: 
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(27) 



which is positive-definite by the convexity of i? as a function of filling. Similarly, for ^ > 2, we choose p'/q' 
p/q+{l + l)/L 

P'\ fp\ T.fN + ^ + '^\ T.f^\ r.f^+'^\ + 1 



ML — - Mi?. 
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(28) 



which is again non-negative by convexity. Thus the intervals of stability of different rational fillings do not overlap. 

It remains to show that the sum of the intervals of stability of all rational fillings completely covers the range 
< /i < m(1)l- According to ([24|) . the chemical potential /i(1)l at which the v = \ state becomes unstable is: 

oo 

m(1)l - ^(n + l)V{ii) - iiVin + 1) 

n=l 
oo 

= 2Y,V{n) (29) 

n=l 
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To compute J^xegif^ni^) ~ I^l{x), we must add together the length of the interval for all possible denominators q, 
multiplied by the number of reduced fractions with this denominator. The function which gives this multiplicity is 
the Totient function (f>{q), which counts the number of integers less than q which are relatively prime to q. Hence the 
combined length of all intervals is 



q—2 n 



Collecting the contributions to V{r) for each r gives 

r 

Using the identity [2^ 



^ </>(g)(r + l)+ 0(q)(r-l)-2^0(q)r 

q\r-\-l q\r—l q\r 



Vir) 



E<^(9) 



(30) 



(31) 



(32) 



gives 



/ = ^ [(r + If + {r- If - 2r2] V{r) 

r 

r 

which is precisely the boundary of stability of the ly — 1 state, establishing the desired result. 



(33) 



IV. AWAY FROM THE CLASSICAL LIMIT: MOTT-HUBBARD TRANSITIONS IN THE STRONG 

COUPLING EXPANSION 

The physics of the Mott transition in one dimension has been studied in the context of classical phase boundaries, 
as the commensurate-incommensurate transition, as well as in quantum mechanics via the Hubbard and extended 
Hubbard models. Basically, as quantum or thermal fluctuations in particles' positions increase, commensurate order 
is destroyed by a condensation of solitonic defects. In the quantum mechanical case, these defects cost potential 
energy, but are favored by kinetic terms, driving a transition at some finite hopping on the lattice. This physical 
mechanism is also responsible for a commensurate- incommensurate transition in systems with infinite-ranged convex 
interactions. However, the fractal structure of the classical ground states leads to a phase portrait considerably more 
complex than in cases studied previously. Here we will examine the structure of this phase portrait. 

The qualitative behavior of our system in the convex regime is reminiscent of the Bose-Hubbard model, with 
commensurate Mott lobes ceding to superfluid states as t increases. It is convenient to treat a state with large q as 
a state with smaller g at a nearby filling in which a crystal of dilute solitons has formed. Hopping tends to liquify 
dilute crystals of solitons: at large separation the inter-soliton repulsion is smaller than the kinetic energy gained 
from delocalization. (The latter grows as 1/r^, and the former as at large separations). The delocalized solitons 
destroy long-ranged spatial order, creating a Luttinger liquid with full translational symmetry. Hence, as t increases, 
the system undergoes a transition from the Mott insulating CGS to a Luttinger liquid state, with larger q states 
liquifying at smaller t. 



A. Strong Coupling Expansion 

To find the position of the phase boundary, we generalize the method of Rcf. and compare the energies of the 
CGS and its adjacent qSS to third order in t using standard time-independent perturbation theory. This approach 
assumes that the phase transition is continuous, so that for a given t > 0, values of fi for which Eqssip-) = EcGsifJ-) 
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constitute the phase boundary. This assumption is well-founded, since the soliton repulsion ensures that the energy 
cost of creating multiple solitons is larger that of a single soliton, thus favoring a second-order transition. 

In a finite system one must account for the repulsion between solitons; here we drop these terms and consider 
the infinite volume limit. In the limit that the solitons are sufficiently well separated that we may neglect their 
interactions, the qSS is highly degenerate and can be expressed in terms of a band of solitonic momentum eigenstates. 
Here we consider only the bottom of the band, which lies at zero momentum. 

To find the transition, we compare the energies (calculated up to third order in t) of the CGS (which is the t — 
ground state) and the qSS. The calculation of the energies is carried out using standard time-independent third order 
perturbation theory. The energy corrections are calculated in terms of the t ^ ground and excited states of the 
CGS and qSS. 

1. Perturbation theory in the CGS 

To calculate these energy corrections, we must consider excitations about the classical CGS generated by pertur- 
bative hopping, and their matrix elements with the unperturbed state. The zeroth order CGS is given by the HUP 
solution \'tlj{CGS)), which is non-degenerate (up to global translations). Hopping creates excited states of the form: 

\4l\x))=b%+,\4'iCGS)) (34) 

At filling ^ there are exactly 2p distinct such hoppings which must be considered: one in each direction for each 
occupied site in the HUP unit cell. 

To calculate the energy corrections, we must calculate matrix elements of the ground and excited states with Hi, 
and the energy differences 5Ei between the corresponding classical ground and excited states. Note that the first 
order wave function contains only terms which can be transformed into the classical ground state by only one hopping 
- in other words, only the excited states of the form and ([57]) have non-zero contributions. The details of the 
perturbative approach are outlined in Sect. IIX A[ here we will only state the results. 

The first and third order corrections to the CGS energy are zero, because the CGS is non-degenerate and hence any 
odd number of hoppings produces a state orthogonal to the ground state. At ly ~ p/q, the second-order correction is 
given by 

^ i=l * 

where AE, = £;f' - E^°^ is the differ ence in potential energies between the ground and the excited state formed by 
hopping from the i*'' occupied site in the ground state configuration. As the ground state is periodic, it suffices to 
calculate these energies for the p distinct particles in the repeated pattern. 

2. Perturbation theory in the qSS 

We now repeat this analysis for incommensurate fillings, which have solitonic ground states. The qSS is given at 
zeroth order by the HUP qSS, consisting of q HUP solitons sufficiently far apart that soliton-soliton interactions can 
be ignored. The energy of the qSS is then simply q times the energy of a single soliton on an infinite lattice. The qSS 
ground states are highly degenerate in this limit, since all translates of each soliton have equal energy. We therefore 
work in the basis in which the hopping term is diagonal, namely: 

l<ss(fc))-E^'''l^^^^(^)) (36) 

X 

where the sum runs over occupied lattice sites, and l^pqssix)) is the state containing one soliton beginning at lattice 
site X. The soliton hops by q sites when a single boson on one of its edges is hopped by one site; hence the position 
X is an integer multiple of qa. The perturbation theory is essentially the same for qSS states containing solitons and 
anti-solitons; functionally the difference between these is in the energy gaps to the local excitations. 

Local excitations, or defects, of the qSS are again generated by hopping one particle away from its preferred position. 
These can be expressed in the form: 

|V.(A:)) =5]e^'=^6(x + r)6t(:, + r±l)|^(°yfc)) (37) 

X 
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The state 4'r{k) describes correlated propagation of a soliton and a defect r lattice sites awav[25f . Note that for each 
r the defect has 2 possible orientations, depending on whether the hopping has been towards or away from the defect. 
Computing the resulting matrix elements gives the energy corrections: 



E 



qSS 



E 



(2) 
qSS 



■'2qt cos(fcqa) 



2qcos{2kqa) 



E. 



(3) 
qSS 



-2qt^ cos{kqa) 



N/q 

cos(2fcga) 



AEr 



cos{2kqa) 



AEr,^_,AEr,^_,+g (A£;,,_J2 

N/2 

-2qt'^ cos{kqa) 



i=l 



1 



AE^.AEr 



(38) 



Here k is the soliton momentum, a is the lattice constant, and q is the denominator of the commensurate filling 
fraction, which appears here because one hopping displaces the soliton by q lattice spacings. The subscripts on Er^^±i 
indicate the distance between particle i and the soliton, and the direction of the hopping relative to the soliton. The 
special distance ri describes an excited state in which an anti-soliton is sandwiched between 2 solitons (for a solitonic 
qSS), or vice versa (for a qSS with anti-solitons) . These states contribute extra terms to the energy corrections of the 
qSS because of the ambiguity as to which soliton is associated with the ground state qSS. 



u/ v„ 




0.005 



0.01 0.015 



0.025 



0.035 



FIG. 1. Perturbative calculation of the Mott to SF phase boundary in the {t, plane, shown here for U = 20. Here the strength 
of all couplings is measured relative to that of the dipolar interaction strength Vc). Each lobe encloses a Mott insulating region 
in which the filling is fixed; the region outside the lobes is a superfiuid of solitons. The 1/3, 1/5, 1/6, 2/7, 2/7 and 3/7 -filled 
lobes are shown here. Every commensurate state of the complete Devil's staircase has a Mott lobe, but the range of hoppings 
over which a state exists falls ofi^ sharply with its denominator. 

Fig. [1] shows the results of the perturbative calculation for selected filling fractions. The t = axis corresponds to 
the classical limit, in which the CGS states comprise a complete devil's staircase: every value of /i corresponds to a 
rationally filled ground state, except for a set of measure 0. The figure shows the resulting Mott lobes: inside each 
lobe the CGS is stable and the system is in a Mott insulating state. The Mott gap vanishes on the boundary of the 
lobe; outside of this region solitons proliferate and the system is in a Luttingcr liquid phase. 
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For any t > 0, only a finite number of insulating states exist; the rest are liquid states with a supcrfluid of condensed 
solitons. The function is no longer a devil's staircase, but rather a piecewise smooth function, with plateaux of 
constant density separated by liquid phases whose density varies continuously with fi. The size of the commensurate 
region decreases sharply with q: states of higher q have both smaller ranges of stability in the classical limit, and 
larger e nerg y corrections relative to the CGS, so that the volume of the corresponding Mott lobe scales approximately 
as l/(7^.[26[ The total volume occupied by liquid states can be estimated from the first-order approximation to the 
Mott lobe boundaries; we find that for small but fixed t the volume of the liquid region scales as approximately 
t^^^. (Details of the calculation arc given in Sect. IIX B[) . At sufficiently large t we expect all insulating states to be 
unstable, and the particle density to vary smoothly with fi. 



B. Bosonization treatment of the phase transitions 



The qualitative nature of the Mott transition can be deduced from existing knowledge of ID commensurate- 
incommensurate phase transitions, which we summarize here. We approach the phase transition from the Luttinger 
liquid side, where bosonization can be used to treat the kinetic term and dipolar interactions exactly. The lattice 
potential can be added perturbatively; when the coupling of the density to the lattice becomes relevant, this signals 
the Mott transition. 

In bosonizcd form, our system is described by a Hamiltonian that is a sum of a Luttinger piece, accounting for the 
kinetic term and dipolar interactions [27| . and a sine-Gordon term which emulates the lattice in the continuum limit. 
We also include a parameter S, which parametrizes the deviation from commensurate filling. Close to filling ^, the 
bosonized Hamiltonian is thus: 



H = — I dx 
2-K 



uK{TTn{x)f + -^(V0(-^))^ +gcos{2q(j){x) - Sx) (39) 



Here </> is related to the density of the physical bosons by 



p{x) 



1 

po V ' 



and H is the momentum conjugate to (f). u and K are the usual Luttinger parameters, expressible in terms of the 
hopping and interaction terms of the original system. The sine-Gordon term adds a periodic lattice potential 

27r 

dxp[x)V cos{ — x) (41) 



Near commensurate densities Po = ~ + the expression (|40p for the density shows that the slowly varying modes 



to the model; g parameterizes the strength of the coupling to this lattice 

Near commensurate dcnsit 
are precisely those for which 

n^q (42) 

where n is the resonant harmonic from the density (1401) . The resulting expression for the coupling (j4ip to the lattice 
is 

V J dx coa{2q(t){x) ~ S) . (43) 

For densities ^ we obtain the same expression, since the state at filling ^ consists of several charge density wave 
instabilities at wave vector q, whose relative positions also become pinned inside the Mott region. 

The Hamiltonian ([39]) is well studied in both the context of the Mott transition [l^ and the Frenkel-Kontorowa 
model of surface interfaces [2[. The upper and lower Mott lobes join in a cusp which is not accurately described by 
the perturbation theory- since 'small' hopping implies that the ratio t/SE of the hopping relative to the energy gap 
to the solitonic states must be small, the range of t over which the pcrturbative treatment is valid decreases with q 
and never encompasses this point. 

The bosonized treatment reveals that crossing the edge of the Mott-Hubbard lobe induces one of two different types 
of phase transitions. At the cusp joining the upper and lower Mott lobes, a constant-density phase transition, with 
(5 = in Eq. of the Kosterlitz-Thouless type occurs. For S the transition is well described by a simple 

two-band model with quasi-particlcs that are gapped in the commensurate phase, and with a density increasing as 
y/p — pc near the transition on the liquid side. 
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Derivation of the Luttinger description near the phase transition 



To derive these results, the Hamiltonian (p9)) must be treated using a perturbation series in g. This is easy to do 
for 5 = 0. One finds that for K < Kc, the cosine term becomes relevant at low energies, and locks the state into a 
crystalline phase commensurate with the lattice. For K > Kc, so-called domain walls (in our case, lattice solitons and 
anti-solitons) proliferate, and the phase becomes a liquid. In the renormalizcd Hamiltonian, this predicts a transition 
from the liquid to the solid phase at < ~ 1/g^ at filling ^, consistent with the perturbative results of Section HV Al At 
the transition, K jumps discontinuously from the universal value Kc = 1/2 (approaching from the liquid state) to 
in the gapped Mott phase. This is the standard signature of a Kosterlitz-Thouless transition. 

A phase transition with (5 = occurs only for the value of ^ for which the model is particle-hole symmetric. For non- 
zero S, the situation is similar, but at constant /i the renormalization of the doping flows to strong coupling before 
the phase transition; hence the transition is not accurately described by perturbation theory about the Luttinger 
Hamiltonian. In this case the transition is best studied using the Luther-Emery solution- that is, rc-scaling the 
problem and exploiting the fact that it can be mapped onto a system of weakly interacting fcrmions (in our case 
actually the lattice solitons) with an upper and lower energy band separated by a gap. 

To do this, we first scale the fields in Eq. ([M)) according to cj) ^ qcp + Sx/2, H = H/q. Then 



H = — I dx 
2tt 



uK{TrIl{x)f + -^(V0(2;) " Sx/2f 



-t-cos(2(/)(a;)) (44) 



where K = Kq^. The C-IC transition occurs at a value of if = or in this picture K = 1/2. Here we will drop all 
~ and work in the scaled system. 

Second, we translate the Luttinger Hamiltonian (j44p back into the language of spinless fermions: 

k 

+92{pR{k)pL{-k) + H.C) + 34 [pRik)pR{-k) + pL{k)pL{-k)] (45) 

This is convenient because the potential cos 2(/) due to commensuration effects has been mapped to the quadratic term 
c^Cl -|- h.c. Thus at K = 1, along the Luther-Emery line where the system is particle-hole symmetric, the interaction 
terms 172 and 34 vanish and we have: 

Hq = VFk{c'i^^.CRk - c\f.CLk) + ^{c^rCl + h.c.) (46) 



This simply describes free massive fermions, of energies e = ±y^v'pk'^ + A^, with the mass gap A is set by the coupling 
of the cosine term. For K ^ 1, Schulz Q showed that a paramctrization can be chosen such that the 54 term vanishes, 
and the coupling of the remaining interaction is proportional to the doping away from commensuration, which vanishes 
at the phase transition. Thus the non-interacting model gives a good description of the physics very near the phase 
transition. 

In practice, the quasi-particles in (|46p may be identified with Hubbard solitons (or distortions of the commensurate 
lattice) which become gapless and proliferate at the transition. In the solid phase, the solitons are gapped, and do 
not occur at sufficiently low temperatures. At the phase boundary the chemical potential crosses the bottom of the 
upper band, and solitons proliferate. Near the transition on the liquid side, this predicts d ^ ^ p — pc with pc the 
chemical potential at the phase transition. 



V. EFFECTS OF THE TRAPPING POTENTIAL 



In practice, any experimental realization of the HUP model will involve a finite sized atom trap, generally of length 
not more than a few hundred lattice sites. In addition, the trapping potential is not generally flat, but rather is well 
approximated as a harmonic potential. Here we will try to address the possible effects of this trap on the system. 

We can estimate the effects of a harmonic trapping potential (present in current cold-atom experiments) using the 
local density approximation - this assumes that the trapping potential is slowly varying enough that it can be simply 
incorporated into the chemical potential, resulting in a spatially-dependent p. Trajectories along the ID chain then 
correspond to cuts in Fig.[T]at fixed t/V[). Thus, at t = we expect different commensurate fillings at different points 
along the trapped chain, with the most stable states (fillings 1/2 and 1/3) occupying the largest regions within the 
trap. However, commensurate states with a period q greater than the lengthscale over which the trapping potential 
varies cannot exist in the trapped system since they violate the local density approximation. The situation is improved 
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slightly when t is nonzero but small, where small islands of Mott states at various fillings p/q are separated by regions 
of superfluid. These fluid regions can interpolate continuously in density between the two commensurate states. We 
will discuss the resulting density profiles for various chemical potentials. In doing so, we ignore the CDW correlations 
of the liquid states, which near filling p/q also occur at length scale q; hence this treatment is a modest improvement 
of the LDA. 



A. t — physics and the LDA 



The simplest approach to determining the density profile in a finite harmonic trap is to use an effective Local 
Density Approximation (LDA) treatment. In other words, given an approximation to the function ^(/i) as calculated 
by @, one can generate a profile of the filling v{x) as a function of position, x, in the trap, treating the trapping 
potential as a spatially varying chemical potential. In this approach, we neglect the fact that each filling p/q requires 
at least q lattice sites to be realized, and that the range of x over which v = p/q may well encompass less than q sites. 
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FIG. 2. Local density approximation (based on the HUP solution for infinite /i) for a trapping potential V{x) = {x — 
Xfyf' / {L / A)a^ — Ho (where L is the total number of sites in the trap) shown for various values of fio- We expect the large 
plateaux of small denominator states to remain as features in the actual finite- volume sohition, but that most smaller plateaux 
will disappear into regions of ambiguous density. 

Figure IV Al shows sample LDA profiles for various values of the trap minimum /xq i taken from a numerical compu- 
tation of i/(/i) cut off at q = 200. For some parameter choices, most of the trap lies in a chemical potential range that 
is strongly locked at a commensurate filling with low denominator, and the LDA gives a plausible rendering of the 
density profile over much of the trap. However the figure also shows clearly that for some regions of the trap the LDA 
is a very poor approximation, as it predicts fillings v = p/q with q much too large to fit into the number of lattice 
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sites over which the state is stable. 

In the classical limit, another approach is to simulate the density profile numerically. We have used simulated 
annealing to generate profiles; Figure W Al shows numerically generated profiles for the same trapping potentials as in 
Figure fVAl Note that in regions of transition between different commensurate states, the exact density is somewhat 
ambiguous, and can vary depending on what technique is used to calculate it. 
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FIG. 3. Density profiles calculated using simulated annealing for the same potentials as in Figure IVAI The multiple lines show 
multiple runs of the simulation, indicating that the states at low-denominator filling fractions are very robust, while states at 
higher values of g (g > 4) are energetically delicate and require simulations over longer time periods than those undertaken 
here. 



Though neither of the techniques discussed above proves particularly illuminating in the transition areas between 
plateaux of relatively small q states, they illustrate some practical features of the HUP system in a trap. First, the 
spatial homogeneity of the density depends strongly on the depth of the trap. Since the range in /i over which the 
lowest q states = 1/2, 1/3...) are stable is much larger than the ranges of stability of higher denominator states, a 
trap with a modest curvature can be 'locked' everywhere at a density of 1/2. As is shifted, however, portions of 
the trap fall outside the range of half-filling, and since all of the nearby states have much larger denominators, their 
ranges of stability will be correspondingly much smaller. In a trap of only a few hundred sites, for the most part 
these states will not be stable over a wide enough range to produce a strong scattering peak. Hence it is the strongly 
locked regions which are the most visible experimentally, and these, fortunately, are well -described both by LDA and 
simulated annealing techniques. 
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B. Spatial Profiles of Atoms in a harmonic trap at finite t 



We have seen that the LDA gives a poor description of the expected profiles for high-denominator filhng fractions 
in a finite trap. However, in practice such states would not in any case exist in any currently attainable experiment. 
Hence in practice we expect density profiles which consist of phase separated regions of density-locked Mott states 
interspersed with compressible fluid states whose density varies continuously with the trapping potential. 

The qualitative properties of the quasiparticlc fluid, at least near the transition, are well described by the Hamil- 
tonian (|46|) . To describe the density of the Luttinger liquid regions, we make a quantitative mapping between the 
2-band model (HH), and the perturbative description of the quasi-particle fluid valid for small t. 

If we ignore the weak interactions between quasi-particles, this effective description has 2 free parameters: the 
effective Fermi velocity in the upper band, and the band gap. By calculating these parameters using the perturbative 
hopping model, a quantitative matching between the Luther-Emery solution described in Sect. II VB H and our system 
can be achieved. We will use this model to generate profiles of particle densities in the trap at finite filling. 

The band gap A is given simply by the difference in energies of the IC states with one extra particle and one extra 
hole. At i = this is just the devil's staircase pattern calculated by @; at finite t it can be extracted from the 
perturbative calculation of Sect. IIV Al 

To match the Fermi velocities, we match the true free soliton Hamiltonian: 

Tjt*) = s\k) [A + fi + t{l- cos{qak))] s{k) (47) 

onto the effective Hamiltonian 

= /(fc)V(«fc)2 + A2s(fc) (48) 

where v is the Fermi velocity in the original 2-band model. 

Both and are quadratic in k near the bottom of the upper band. The effective model describing the 
dynamics of particles in the upper band only is obtained in both cases by linearizing the spectrum about kp in the 
upper band. This gives: 

HW = s''{k) [A + /i + t{l - cosiaqkp)) + kaqtsinlaqkp)] s{k) 

Hi^) = v^kFk/^{vkFf + A2 (49) 

Of course /.i is such that the constant terms in the first equation cancel, and both energies are linear in k. Hence we 
can match: 



v'^kp/^/ivkFY + A2 = qatsmiqaKp) (50) 

To lowest order in kp (valid at low quasiparticle densities) this gives: 

{qafAt (51) 

Wc may use this result to calculate the density of solitons in the IC region near the phase transition. In particular, 
if we ignore quasiparticle interactions, we have 

2nvJ^ Ve2 - A2 
= A.^^I?^aJ (52) 

where in the non-interacting model A = is the chemical potential at the transition, which can be estimated using 
the perturbative calculation of the phase diagram. (Interactions will renormalize the gap in principle, though in our 
case we calculate the gap perturbatively, and this result should be accurate for sufficiently small t.) 

Substituting in the value of v obtained in equation (|5ip , and accounting for the fact that each soliton has a charge 
density of | relative to the background charge density of the lattice, we obtain the charge density near the C-IC 
transition: 

p{q) = -J-—^/J]^-J4 (53) 



This gives the desired expression for the density of quasi-particles as a function of the parameters of the original HUP 
Hamiltonian. 
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FIG. 4. Density profiles in a harmonic trap at tjV = 0.001. Stable plateaux can be seen at 1/2,2/5, 1/3, and 1/4 filling, 
separated by Luttinger liquid regions of continuously varying density. 



C. Density profiles 



Wc now use the results of the previous section to estimate the density profiles in a realistic trap, by using a finite 
t LDA. This scheme first identifies the commensurate regions of the trap; since for the values of t considered here 
all but a few very low-denominator fractions are unstable, the problem of patterns too short to fit into the allotted 
number of lattice sites does not arise. At the edges of each commensurate region we use ([53|) to predict the local 
density profile; a polynomial interpolation is used to join the various liquid regions. 

Figure IVC] shows typical charge density profiles in a trap ijV ~ 0.001. Commensurate plateaux at filling fractions 
1/2,2/5,1/3, and 1/4 can be stabilized, depending on the chemical potential at the bottom of the trap. Between 
these plateaux we see regions of Luttinger liquid. 

Figure lV CI shows density profiles for tjV ~ 0.02, the value in principle attainable by experiments on polar molecules. 

Though we do not claim to predict accurately the density profile deep in the liquid region, where both soliton 
interactions and lattice commensuration effects have a strong impact on the solution. Figure IV CI gives an accurate 
representation of the commensurate regions and their immediate vicinity. Experimentally speaking, the most inter- 
esting features of these profiles are the spatially separated commensurate plateaux, which even at tjV = 0.02 can 
cover a significant fraction of the trap's volume, and hence should be experimentally visible in the structure factor. 
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FIG. 5. Density profiles in a harmonic trap at the experimentally realizable value tjV = 0.02. Here only the 1/2 and 1/3 filled 
commensurate states are stable. 



VI. DEPARTURES FROM CONVEXITY 



The PUH CGS are the classical ground states so long as the potential is everywhere convex. Since the on-site 
potential U is tunable experimentally, it is interesting to ask what happens to these states as \J is lowered away from 
convexity and double occupancies begin to form. Of course, as U is lowered still further, triple and higher occupancies 
will also form. However, as the barrier to triple occupancies is 3C/, instabilities towards triple occupancy at a given 
filling will set in at approximately one third the value of U for instabilities to double occupancy. For U in the range of 
relevant for interesting physics about the 1 /2-filled state which we will describe below, for example, triple occupancies 
will not be favored at any filling fraction. Hence we will not analyze triply occupied states here — the doubly occupied 
regime has enough challenges of its own. 

As is deformed away from convexity, a series of thresholds exists, at values of C/q decreasing monotonically with 
the density of bosons. At each (commensurate) density, we consider two thresholds. The lower threshold is where 
double occupancies start to form spontaneously in the CGS, and a new non-HUP classical ground state takes over at 
this commensurate filling. While our primary interest is in ground state transitions, much insight is gained by also 
computing a second, upper threshold at a given filling. Past this threshold, U = Uc''^^\ a particle added to the state 
goes in as a double occupancy instead of fractionalizing into q solitons. We will see that these two thresholds allow 
us to understand many striking features of the phase diagram. 

(red) for states with q < 15 in the vicinity of half- filling. The values of 



Fig. m plots [/i"?^^^ (blue) and C/< 

(qSS) 



(CGS) 



shown there are obtained by numerical minimization in the sector with one added charge at the specified 



filling. The values of U, 



(CGS) 



are obtained by numerical minimization over configurations at the specified filling that 



20 



contain exactly one double occupancy. The latter is the correct answer for all v = p/q for which p and q are not 
both odd. In such cases the double occupancy and its surrounding charge rearrangements give rise to even moments 
starting with the equivalent of a quadrupole. Consequently, double occupancies repel at all distances and enter via 
a continuous transition at the computed threshold. However, when p and q are both odd, double occupancies have 
a dipole-like moment, causing them to attract at long distances. The transition in this case is first order, and the 
true f/i'^'^'^' lies above our numerically determined value. We ignore this gap due to first-order effects here, as we 
do not expect it to be very large. Indeed, a relatively straightforward calculation shows that the minimum in the 
dipolar potential between two such defects (see Sect. IIX C[) places the two defects at least a distance q apart; hence 
this energy gap decreases at least as This is a small perturbation in states of filling p/q where p and q are both 
large. 

This section outlines the energetic arguments for the locations of these thresholds and the nature of the new ground 
states, and discusses the interesting features of Fig. [6l We will end by exploring what can be said about the phases 
deep in the non-convex regions - which will lead to an interesting new series of states which will form the subject of 
the last two sections of this paper. 
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FIG. 6. Numerically calculated values of Uc (red) and Uc (blue), for a selection of filling fractions v. Quoted values 
of U are measured relative to Vq. A segment of the y-axis between ui'^'^^\l/2) and L^i'^'^' (1/2) has been removed for better 
resolution of the rest of the phase diagram. The green and yellow boxes, bordered by black dots, indicate the approximate 
regions where the ground states consist of double occupancies in the 1/3 and 1/2-filled states respectively. Between the shaded 
regions the ground states consist of double occupancies on other, higher-denominator states. 



A. Effective convexity and thresholds for double-occupancy formation 

We begin by understanding where each CGS becomes unstable to forming double occupancies. Let ujf^'^^\v) be 
the thresholds at which the PUH ground states give way to ones with at least one double occupancy — these are marked 
as the red points in Fig. ([5]). Observe that these thresholds increase monotonically with v, Uc'^'^^\i^') > ui'^'^^\v) 
for v' > V. 

To understand this monotonicity, along with the approximate locations of these thresholds, we considering the 
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implications of convexity at a given filling. Sufficient conditions for convexity [l| are that, for all x, 

i(F(0) + y(2x))>y(x) or U>fvo(^-^^ (54) 

For convexity to hold everywhere, (|54p must hold for x = 1; below this threshold double occupancies may occur. 
However, when perturbing about a given convex solution at fixed fi, solutions will be stable approximately until U 
violates (|54p for a; = r,„, the minimal inter-particle distance. This implies that states with lower filling fractions are 
more stable against double occupancies, as the potential gain in lattice energy from doubly occupying a site is smaller. 
This is in contrast to the stability of the commensurate states as t increases, where the denominator of the filling 
fraction determines stability. 

For example, consider the states v < 1/2. In reality, these states contain no pairs of particles separated by a; = 1, 
so the convex solutions should be stable until Uq violates ((54|) for x = 2 - i.e. for Uq > .235. More generally, we 
expect the Hubbard solutions to remain the ground states so long as Uq satisfies Eq. (|54p for x > r^, the minimum 
inter-particle spacing. 



Arguments for the approximate locations of the transitions 



For fillings v = 1/q we can understand this more rigorously by considering the energetics of nearest neighbors. (The 
potential 1/r^ falls off quickly enough that nearest-neighbor interactions dominate the energetics in this case). Since 
all particles in the initial CGS have the same separation, forming a double occupancy in the CGS results in replacing 



{q + l)V{q)^{q + 2)V{q + l) + U^ 



(55) 



Removing a particle entirely changes the potential energy by {q + l)V{q) — qV{q + l). The double occupancy, however, 
adds 2 NN distances of g + 1 to the lattice: one of these was originally a NN separation of q; the other was a NNN 
separation of 2q - which wc omit here as it is a higher-order term. 

Equation ([55]) shows that the scale at which double occupancies become energetically favorable in the CGS is 
determined by q: if ( I54p holds for x — q, then double occupancies will not form. Indeed, writing 



q + l 



(g + l) + (l 



1 



)0 



(56) 



we see that for V convex. 



Viq) = V{^{q + 1) + 0) < -^V{q + 1) + (1 - ^)^^(O) 



9+1 

or {q + l)V{q)<qV{q+l) + UQ 



(57) 



Hence if V is convex at the length scale set by q, double occupancies cannot be present in the classical ground state. 

Similarly, we can gain an intuitive grasp for the energetics of ui'^^^'^ by considering only nearest-neighbor interac- 
tions, and repeating the above analysis. If a soliton is formed in the qSS, the q solitons have replaced 

{q-l)V{q)^qV{q-l) (58) 



whereas inserting a particle as a double occupancy induces 2 nearest neighbor distances of q. The net NN difference 
in lattice energy between adding a particle as solitons and adding a particle as a DO is thus 



5E = (q + l)V{q) + Uo- qV{q - 1) 

In this case, DO's will not form if Eq. is satisfied for x = q — 1. 

Thus, we expect that the CGS and qSS will become unstable to double occupancies at approximately 



(59) 



U <Uc{CGS) 
U < Uc{qSS) f 



15^0 1 

15^0 1 

8 [q-lf 



(60) 



In short, states of higher density will become unstable to double occupancies at larger values of J7, and at a given 
filling there is a finite gap between ui''^^^ and u!f^'^^\ 
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What about more general filling fractions? Repeating the above analysis, we see that for fillings p/g, p > 1, nearest 
neighbor distances are r, r + 1. Double occupancies in the CGS replace (r + l)V{r) by (r + 2)V{r + 1). Likewise 
the energy cost of a soliton is (r + l)V{r) — rV{r + 1), while inserting an extra particle as a double occupancy adds 
2V{r + 1)-^^ Hence the nearest neighbor energy cost of forming a double occupancy in the CGS and qSS is: 

CGSSE^ (r + 2)V{r + 1) - (r + l)V{r) + Ua 

qSSSE^ {r + 2)V{r + 1) - (r + l)V{r) + Uq (61) 

Hence in general, we expect the approximate value of C/q at which these transitions occur to be set by the smallest 
nearest-neighbor distance in the CGS. However, the scale of the difference between Uc'^^^^ and Uc^'^^^ is determined 
by further neighbor interactions - in fact, not until p*^ neighbors are included in the energetic calculations is the set 
of allowed separations different in the qSS and CGS. Hence this splitting will be much smaller in this case than for 
the \/q -filled state with q = r in Eq. ([6T|l . 

We now begin to understand the basic features of Figure [6] To a first approximation (obtained by considering only 
NN distances), we have: 

C/f ^(i) = [/f (62) 

and, for all v with ^ < v < 

U?^'{.)=Uf'[.) = Uf'[^) (63) 

That is, the value of Uq at which double occupancies first appear is approximately decreasing with the filling. There 
is a large gap between its values for CGS and qSS states at filling fractions i; the jump for fractions > 1 is much 
smaller. 



2. Form of the doubly occupied ground states and computing the exact loci of transitions 

To calculate the exact loci of the transitions requires finding the exact occupation pattern in the doubly occupied 
states, and computing its energy relative to that of states without double occupancy. The second part of this task 
is easy to do numerically, if not analytically; here we discuss the energetically optimal configurations with double 
occupancy in both qSS and CGS states. 

qSS states: 

We begin with the somewhat simpler qSS states. To find the optimal configuration, we first search for the optimal 
position in the repeating CGS pattern at which to add a double occupancy. Intuitively, this will be at the occupied site 
with the lowest local charge density. Second, we ask whether adding a charge at this site creates further distortions of 
the CGS state. Though we will not prove that this is the globally optimal configuration, the resulting configuration 
will give the least distortions of the CGS configuration for a given local charge density, and hence should be the 
ground state. 

First, how do we identify the locus of lowest local charge density? It is useful to consider a few simple examples: 

• For states of the form 1/q, there is only one site per unit cell - hence no freedom in where to place the double 
occupancy. Using the notation of Section FlH A 21 we see that states of filling 2/q have occupancy patterns 
... rr + lrr + 1..., and again both sites in the unit cell have the same local charge density. 

• At filling 3/g, the occupancy pattern is either ... r r r + 1 r r r + 1..., or ... r r + 1 r + 1 r r + 1 r + 1.... In the 
former case, the DO cannot sit between the two NN distances of r; in the latter, it must sit between the two 
NN distances of r + 1 . 

• In the state 5/13 = ...2 3 2 3 3..., the DO must sit between the two NN distances of 3 - which we will henceforth 
denote as the 5*'' site in the unit cell, using the convention that wc count from left of the first interval shown. 

• In the state 5/12 = ...2 3 2 2 3..., the DO must sit on the first or second site in the unit cell (these being 
equivalent, up to refiection of the entire pattern). In this case, considering only nearest neighbors would suggest 
that it could sit anywhere but the 4*'' site. However, occupying sites 3 and 5 leads to second neighbor distances 
of 4, 5, instead of 5, 5 for sites 1 and 2, making these energetically preferable. 
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This suggests the foUowing simple algorithm: begin by considering the sum of the two nearest-neighbor distances 
Si(i) = rn + r to the left and right of each site in the pattern. If one of the Si(i) is less than the rest, place the DO 
here. If not, generate the set of sums of right and left second neighbor distances: S2{i) — rij + Vi^ij + ri^r + n+i,r- 
If this set has a unique largest element at site z, place the DO there. If not, proceed to the set of sums of right and 
left third neighbor distances, and so on. One can proceed in this way up to p— 1''* neighbor distances; sites for which 
Sfe(i) = Sfe(j)^~i[ have identical local charge density - as in the example above at = 5/12- and hence are related by 
reflections and translations of the unit cell. 

Having found the optimal site at which to add the DO, we must then ask whether inserting charge here results in a 
further re-arrangement of the background charge of the lattice. Basically, dipolar 'charges' on nearby sites are repelled 
by the extra charge at the DO. Charge will be forced away from the doubly occupied site when the potential energy 
gain in doing so is greater than the cost of compressing the state by the corresponding amount in the remainder of 
the lattice. When such compression is favorable, it will result in particle-like solitons forming near the DO and being 
repelled to infinite distance. In practice, this means that displacing one charge, say at the right of the DO, one lattice 
site outwards will push all charges to its right one site outwards as well. This creates a single soliton at oo, as we 
have lengthened one p^^ neighbor distance from g to q -I- 1. 

It is straightforward to calculate the relevant energies. The lattice energy associated with a single soliton given by 
Eq. (IH: 

^Ei^l/q J2 [{rp + l)V{rp)-rpV{rp + l)]+J2mVinq-l)^{nq-l)Vinq) + ... . 
P7i0(mod q) " 

We must compare this to the energy gained by creating an extra hole between the DO, and site i, a distance d to its 
right (say). This energy has two contributions: one infinite sum for the change in interaction of site i and all particles 
to its right with the double occupancy, and one double infinite sum for the change in interaction of site i and all 
particles to its right with all occupied sites to the left of the DO. This gives: 

oo ^ 

Ai^2 = E 



(xk ~ Xi + d + 1)3 [xk ~ Xi + d)3 



. {xk + Xj - Xi + d + 1)3 {xk + Xj ~ Xi + dY 

j — ~oo k—i ^ / \ J 



(64) 



where Xk here denotes the position of the fc*'' particle on the chain, relative to an arbitrary origin. The sum over j 
here includes the doubly occupied site; the first line of Eq. (IMl) counts only the extra interaction due to the second 
particle occupying this site. 

The important point here is that the sum of the second line of Eq. and the soliton energy is positive-definite 
by convexity of l/r^. Further, their sum depends only on which particle in the unit cell lies at site i, rather than 
the value of d directly. The gain in energy from moving away from the DO (first line of Eq. conversely, falls of 

rapidly with d. Hence if d is sufficiently large, the repulsion of the DO will be too weak to push the charge outwards. 
Further, once we have found a site for which A£^i -I- AE2 > 0, no charge further from the DO than this site will 
be displaced. It is thus straightforward to compute AEi + AE2 for the sites close to the DO to determine which 
charges will be pushed outwards. The algorithm begins with sites closest to the DO; if AEi + AE2 < 0, the charge in 
question and all charges further from the DO are pushed outwards by one site, creating a single soliton infinitely far 
away. One then moves outwards to the next closest charges, and repeats the process until a site is reached for which 
AEi + AE2 > 0. 

CGS States 

Forming double occupancies in the CGS is qualitatively different, since we think of re-arranging the existing dipoles, 
rather than adding a new one. In this case it is simplest to think of forming the DO by pinching charge inwards 
towards an unoccupied site on which the DO will form. 

The pinching operation involves picking a site Xd on which the DO will be formed, and moving some number of 
particles (initially at positions {Ri}) to the left of Xd rightwards, and some number of particles to the right of Xd 
(initially at positions {Li}) leftwards. 

For example, the 1/3 state introduces {Li} ~ {!}, and {Ri} ~ {2, 5}, relative to Xd- A DO can be formed in the 
1/3-filled state as follows: 

i:1001O01001001 

Pinch Li,i?i: 1000200001001 

Pinch i?2:1000200010001 (65) 
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Similarly, for the 2/5 state, {L,} = {1,4}, and {i?,} = {1,4} 

-:101001O100101001 

5 

Pinch Li,i?i:1010002000101001 

Pinch i?2:1001002001001001 (66) 

Hence again, we must first find the optimal site onto which to pinch the charge, and then ask whether other charges 
will move after the DO has formed. The lowest energy final charge configurations are as symmetric as possible: if the 
local charge density is lower on one side of the DO, this inhomogeneous distribution will exert a force on the excess 
charge of the doubly occupied site. Thus we find the optimal site by searching for holes which sit at symmetry centers 
of the configuration. 

This approach reveals an important even-odd type effect in forming DO in the CGS. If g — p is odd, every hole 
in the unit cell has an equal number of holes in the pattern to its right and to its left, and we will show that it is 
possible to choose a hole about which occupied sites are distributed symmetrically. Pinching in these states will not 
alter this symmetry, and hence two DO's repel. If g — p is even, however, such a site does not exist and the final 
charge distribution is not symmetric. This leads to a dipolc-like moment for each double occupancy, causing these to 
attract at long distances. 

As a simple example of how this works, consider the states at 2/5 and 3/5 filling. The occupancy patterns are: 

i/ = 2/5 : ...1010010100... 
= 3/5 : ...1011010110... 

(67) 

The important difference between the two is that at i' ~ 2/5, there is an odd number of holes in the unit cell - and 
hence somewhere in the unit cell there is a pair of occupied sites separated by an odd number of holes. The distribution 
of charge required by convexity is such that the hole at the center of this interval is a center of symmetry of the charge 
distribution. The double occupancy can then be formed by pinching particles onto this center of symmetry. The 
ensuing charge distribution is therefore symmetric. At v = 3/5, the center of symmetry of the pattern occurs on an 
occupied site, and hence there is no way to pinch the charge completely symmetrically. The result is a configuration 
which has a dipole-like moment due to this charge asymmetry. 

Indeed, in general HUP solutions will have a center of symmetry, since this is the most homogeneous way to 
distribute charge. For general p/q with p — q odd, the unit cell contains an odd number of holes, and this center of 
symmetry consequently lies on an unoccupied site. For p ~ q even, conversely, it lies on an occupied site. Thus the 
HUP solutions also result in the even-odd effect observed in the phase transitions to doubly occupied states. 

This pinching construction gives a simple algorithm for constructing DO in the CGS: first we find a hole which sits 
at a center of symmetry, and pinch the two nearest bosons onto this site. This can only create attractive interactions 
from the remaining charges towards the DO. We then iterate through the neighboring particles and compute the energy 
of moving each one inwards. (Here one must bear in mind that for less dense CGS states, it may be energetically 
favorable for particles to move inwards by several lattice sites). If this energy is negative, we move the particles 
inwards and proceed to the next nearest occupied sites. If the energy is positive, no further re-distribution of charge 
will occur. 

Hence for both CGS and qSS states, a simple algorithm exists to find the optimal distribution of occupied sites, 
given that a double occupancy will form. This allows us to compute the exact difference in potential energy between 
configurations with and without a single double occupancy - and hence calculate the threshold values of Uq at which 
double occupancies begin to form. The result is the classical phase portrait Fig. [HI 



B. Arguments for monotonicity of Uc with filling fraction 

From our above considerations, at fixed filling, the CGS is always more stable against forming double occupancies 
than the qSS. Indeed, a simple estimate suggests that the gap separating these instabilities is approximately the 
difference in energies for adding and removing a particle to the state- which is equal to the interval in fj. over which 
the state is stable at infinite U. We now use this to understand an intriguing property of Fig. [5] - namely, the 

(CGS) 

monotonicity of Uc with i^. 

The monotonicity can be understood as a simple consequence of the fact that ui'^^^\v) > ui^'^^\v). Consider a 

filling v' > v. As v' can be constructed by adding charge to v, we conclude that m ' {v') > ' [v); the latter 

corresponds to the threshold at which adding an extra filling v' ~v in the form of the solitons of v loses out to adding 
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it in the form of double occupancies on top of the PUH state at v. As filhng factors only slightly greater than v 
involve a dilute addition of charges, we further conclude that \i-aiyi^^+ u!f''^^\v') — Uc''^^\v). 
A somewhat more involved argument shows that lim^/_>^- ui'~^'^^\v') = Uc^^^\iy) . We have: 

hm [/i''^^)(:.')=C^f''^^M (68) 

since the state at filling v ~ S is the state at filling i/ with a density 6 of hole-like solitons. As this density decreases a 
single particle added as a double occupancy will induce a charge configuration increasingly similar to that of forming 
a double occupancy in the p/g- filled state (which can be thought of as forming q hole-like solitons by removing a 
single particle, then re-inserting this particle as a double occupancy and letting the charge settle into its optimal 
distribution.) Further, the gap between ui'^^^\v') and U^(j'^^\v') vanishes as the denominator of the state v' goes 

(CGS) 

to infinity, implying continuity of Uq from the left. 

In short, our considerations so far imply intricate behavior for the location of the initial ground state instability, 
namely that Uc^^^\v) is a monotone increasing function of v on the set of rationals with a discontinuity at each 
rational value of i>: 

hm {7fG^)(:.') >t^f''''ni^) = lim C/f^^)(^')- 

The first relation is a strict inequality, as the scale is set by the gap between ulf"^^\v) and Uc'^^^\v)^ which is finite. 
The second is equality, because the gap is set by lim^^o ui''^^\v + 6) — ui^^^\v + 5), which is 0. 



C. Structure of the doubly-occupied regime 

The previous sections have outlined some rather stringent constraints on the thresholds ui'^^^'^ and Uc'~^'^^\ which 
give significant insights into the nature of the phase portrait of the classical system as Uq is decreased: monotonicity 
locates the transitions of higher denominator states relative to those of lower denominator states. Thus we can 
understand the coarse features of this instability by considering first the most stable states (of denominator q < qo, 
for some go small enough to allow the exact thresholds to be computed), and deducing the expected behavior at 
fillings close to these. 

We now turn to the question of what can be said about the phases with U < Uc- At any given filling, tracking 
the evolution of the ground state with decreasing Uq after double occupancies have been introduced is a problem of 
considerable complexity. Here we use the ideas developed thus far to identify a family of regions in the {iy,Ua) plane 
where simpler descriptions emerge — these are indicated, in two simple cases, by the shaded regions on the figure. The 
basic idea is that once the v = p/q qSS becomes unstable to double occupancy, any particles added to the v ^ p/q 
state will be added as double occupancies, since these repel less strongly than solitons (see below). Hence at first sight 
we expect, in the region ui'^^^\p/ q) > U > Uc'^^^\p/q), states of filling v > p/q to consist of double occupancies in 
the V = p/q state. At rational fillings the double occupancies will arrange themselves in a crystal thus generating a 
commensuration distinct from that of the underlying p/q state — we will refer to these as doubly-commensurate states. 

A detailed discussion of why, for Uo < U^'^^^\i^), increasing the filling fraction can only form new double occupan- 
cies, and never new solitons, is presented in Sect. IIXDI It is useful, however, to understand the simple physical reason 
underlying this: adding particles as solitons results in a charge distribution that is maximally spread out in space; 
inserting them as double occupancies produces a maximally localized charge distribution. If the change in density is 
infinitesimal, the solitons or double occupancies will be infinitely far apart, and thus by definition if Uq < ui'^^^\v) 
these extra particles will enter the ground state as double occupancies. If the change in density is finite, we must add 
some number of particles per unit length on the lattice. Since the potential l/r^ falls of rapidly in space, the repulsion 
between two added charges in a given distance is smallest when the extra charges are as localized as possible -that 
is, when both form double occupancies. Hence as the density increases, so does the energetic payoff of forming DO, 
rather than solitons. 

Our discussion so far suggests that for all Sv > 0, states at filling v + Sv with ui'~^'^^\i') < Uq < u'f^^\v) consist of 
the doubly commensurate states described above. However, as the density of added charge increases, the parent state 
itself becomes less stable to forming extra double occupancies. This can lead to transitions in which the structure of 
the CGS collapses to a crystal of double occupancies over a background of significantly smaller filling. 

We have carried out a simple analysis of the location of this instability for the 1/2 and 1/3 plus double occupancy 
regions in Fig. ^ at selected fillings and these are marked by the black dots in the figure. We do not know of a 
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general algorithm to compute this threshold at arbitrary filling, but the principal is easy to illustrate for the 1/2-filled 
state. Consider the following two states at filling 2/3 = 1/2 + 1/6: 

... 2 1 1 2 1 1 ... 

...2 2 2 2 0... (69) 

The first state will be stable for Uq < Uc'^^^\l/2), but sufficiently large. However, the energy gained by pinching 
this state to form the second state is greater than the energy gained by simply forming a DO in the 1/2-filled state, 
as the two particles that are pinched are moving away from sites with charge 2, rather than charge 1. The black dots 
in the figure are calculated by calculating the difference in energies of such simple configurations. 

It is interesting to notice that the black dots do not behave monotonically with filling fraction. To understand why, 
consider the following two states at filling 3/4: 

... 2 1 2 1 2 1 ... 
... 2002002020020020... (70) 

Because of the greater density in the first state, the best we can do by pinching is ensure that 2/3 of the double 
occupancies sit at least 3 lattice sites apart. Hence once at least half the sites are doubly occupied, the potential gain 
due to forming extra double occupancies is actually less than for smaller fillings. 

In summary. Fig. |5] gives an accurate picture of the locations of the initial phase transition due to deforming Uq 
away from convexity for arbitrary fillings. We have constructed the general form of the resulting states with double 
occupancy. As Uq is decreased even further, our analysis suggests that there are further transitions to states with a 
higher density of double occupancies, whose structure we have not systematically understood. While we have found 
sizeable regions which can be described as simple descendants of the 1/2 and 1/3 states, we are not at present able 
to estimate the sizes of analogous regions for higher denominator fractions. Of course, to have a full solution of this 

( CGS) 

would be equivalent to tracking the evolution of each as C/ is decreased from Uc 



VII. INTERESTING PHENOMENA IN THE NON-CONVEX REGIME 

In Section IVI[ we have expended considerable effort understanding the various phase boundaries in the classical 
system as Uq is decreased away from the convex limit. We now focus on a particularly interesting region of this 
phase diagram - namely, that in which the ground states are derived by adding double occupancies to a parent CGS 
configuration. We will first discuss the classical limit of these phases, showing that these contain a re-scaled version of 
the devil's staircase. Focusing on the devil's staircase near 1/2 filling, we then consider the effect of adding hopping 
to the mix, and find a phase diagram with super-solid like regions in which a Luttinger liquid of double occupancies 
co-exists with a commensurate 1/2-filled background. 



A. A new staircase 



The discussion of Sect. IVII has led us to the doubly-commensurate states in the (ly, Uq) phase diagram: we remind 
the reader that such states are constructed by periodically doubly occupying some fraction of the sites in a CGS. We 
now discuss the ground state configurations of these double occupancies, and show that in at least some cases these 
doubly commensurate states can form a devil's staircase of their own. 

To explore such states, we will work in the parameter range where increasing density effectively adds double 
occupancies to a parent CGS state, without forming DO in the parent state itself. We have argued in Sect. I VI CI that 
at least for the 1/2 and 1/3 filled states, this description is apt over sizeable regions of the phase diagram. Though 
we have not calculated the lower thresholds for higher denominator states, we expect that a similar description holds 
for arbitrary parent fillings ly and Uc'^'^'^^ < Uq < Uc'^^''^\ at least over a modest range of densities. 

First consider states constructed from double occupancies in the 1/2-filled state, which exist in the region shaded in 
yellow in Fig. ([5]). The energetics of such states can be divided into a) the constant interaction of the parent 1/2-filled 
PUH configuration with itself, b) the constant interaction of the added charges, irrespective of their location, with the 
parent 1/2 filled configuration and c) the interaction of the added charges with themselves. This last part involves an 
interaction between the added charges which is convex again and thus leads to PUH configurations sitting on a lattice 
with a doubled lattice constant. The energy cost of adding a single double occupancy is U + Vd, Hence the doubly- 
occupied sites comprise a Devil's staircase with jj, fi + U + Vd, and the widths of all intervals decreased by a factor 
of 8. Here = | '^^ interaction energy of each double occupancy with the underlying 1/2-filled state. 
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At fixed U, this staircase is complete over the range of fillings for which increasing the particle density infinitesimally 
does not induce 'excess' double occupancies to form in the half-filled background lattice. In the case of the 1/2-filled 
state, for U sufficiently close to the upper cutoff this gives a complete staircase on 1/2 < < 1. 

Similar structures exist for all l/g-fiUed states in the appropriate range of U. As mentioned before, we do not, at 
present, understand the situation for doubly commensurate descendants of general rational fillings. 



B. Supersolids 

Thus far our considerations away from the convex limit have been purely classical. But we can equally consider 
states obtained from these modified classical states upon the introduction of hopping. Specifically, let us consider the 
fate of the doubly commensurate descendants of the PUH 1 /q states considered above. 

In a manner entirely analogous to the problem with which we began this paper, the superlatticc of added charges can 
melt via the motion of its solitons as t is increased resulting in a phase transition between the doubly commensurate 
state and a "super-solid" like phase in which the background l/g-fiUed CGS coexists with a Luttinger liquid. This is, 
in a sense the d = 1 version of the supersolid in higher dimensions, but it is worth noting that the d = I version in 
our problem exhibits a more divergent CDW susceptibility than superfluid susceptibility as T — 0. 

To get a more quantitative account of these new phases, we may repeat the strong coupling treatment above. 
Fig. ([7]) shows the phase portrait at intermediate values of U near i/ = 1/2. The black line traces the infinite- [/ Mott 
lobe, over which the background 1/2-filled state is stable against forming solitons. The red line shows the threshold 
at which it is energetically favorable to add a single double occupancy to the 1/2 filled state. The blue curves show 
the positions of the Mott lobes for the doubly commensurate states. The presence of double occupancies stabilizes 
the 1/2-filled state against proliferation of solitons, so that the background remains commensurate at least within the 
infinite-[/ 1/2-filled Mott lobe, shown in black. This 1/2-filled super-solid phase has also been shown to exist in an 
extended Bose-Hubbard model with second- neighbor repulsion [30|; these numerical results are consistent with the 
phase portrait shown here for small t/V. 



1. Perturbation theory m the super-solid regime 

The phase diagram in Figure[7]is obtained using 4*'* order perturbation theory in t. As the calculation of the energy 
corrections is rather involved, we have included the details in Sect. IIX A II Here we will discuss the form of these 
corrections, and the qualitative features of the phase diagram. 

The calculation is similar to that described in the convex case in Sect IIVI with one key difference: all ground states 
must have holes at every second site. Adding a charge to a doubly commensurate crystal produces an extra double 
occupancy - which in turn produces solitons in the doubly occupied sites. That is, as we do not allow triple occupancies 
to form, and it is energetically unfavorable to place extra particles on the unoccupied sites, we may only add particles 
by forming solitons in the doubly-occupied superstructure. At filling 1/2 + p/2q, for example, the relevant ground 
state is given by the qSS of the p/q filled state, with all distances stretched by a factor of 2. Hence only even-order 
corrections to the energy are no n- vanishing: an odd number of hoppings produces a state with at least one particle 
on a site that is vacant in the parent half-filled state. 

As discussed in the Sect. IIX[ we can therefore describe the corrections in terms of an 'effective' hopping Hamiltonian, 
H'^'^'^ , which describes the hopping between occupied sites. H'^'^'^ contains both a 'trivial' correction, in which a boson 
hops off a given site and back on, and a 'hopping' term in which a boson hops from one occupied site to another. The 
second-order correction to the energy in this case is: 

E^^^ =-{i'o\Hff\,Po) (71) 

The trivial part of Hi gives approximately the same correction for commensurate and solitonized states (except 
for small differences in the energy denominators). The hopping term, on the other hand, has a first-order matrix 
element with solitonized states, but not with the CGS. Thus the situation is analogous to that of qSS solitons in the 
convex limit: solitonized states have a negative energy correction at leading order in perturbation theory which their 
commensurate brethren do not. Thus as t increases the double occupancies undergo a commensurate-incommensurate 
transition, as seen in Fig. [T] 
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FIG. 7. Phase portrait in the vicinity of doubly commensurate states about 1/2 filhng at U = Vo- Again, fi,t, and U in the 
figure are measured relative to Vo- The black lines indicate the boundaries of the Mott lobe at U — 20. The red curve shows 
the chemical potential at which it becomes energetically favorable to add particles to the half-filled state as double occupancies, 
for U = Vo- The blue curves show some Mott lobes of the doubly-occupied staircase region; the region between these and the 
black lines (which delineate the region of stability of the 1/2-filled state at infinite U) is a super-solid state. 



VIII. CONCLUDING REMARKS 



The infinite range of the dipolar interaction does produce, as promised, intricate phase diagrams for the one 
dimensional dipolar bosonic gas. Particularly striking are the singular staircase functions that showed up in our 
analysis in three different settings: in the v{^) curve for the exactly dipolar classical problem, in the function U^^ 
which marks the instability of the PUH states when the onsite U is tuned down and in the v{ijl) curves in selected 
regions of the [v, U) plane. The other main set of results pertain to the presence of a large, indeed, infinite number of 
transitions between Mott crystals and Luttingcr liquids or supersolids. The challenge of observing some of this physics 
in cold atomic gases is not trivial — the major obstacles are getting a reasonable simulacrum of a one dimensional gas 
of infinite extent. On the positive side, the control parameters we study here are eminently tunable. 



IX. SUPPLEMENTARY MATERIAL 



A. General Perturbation Theory 

Here we review the perturbation theory pertinent to treating the hopping terms in the general Hamiltonian, and 
in particular hopping in the super-solid state. 

Let Hq be the unperturbed Hamiltonian, with cigenstatcs of energy Ef'\ The perturbation 5Hi induces 

corrections to the wave function and energies; the i*'' order corrections are denoted jV']*''), i?]*^. When the lower 
indices are omitted, we will be referring to corrections to the ground state energy and wave-function (which is all 
we're interested in here.) 
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Perturbation theory can be summarized by the following equations: 

1^) = 1^(0)) + <5|^(l)) + ^2|^(2)) + 53|^,(3)^ ^ 
E = + + J2i5(2) + ^3^(3) ^ 

{Ho + SH,)\i^)=E[^j) (72) 
We stipulate that (V'o^'' I ''/''■'■') = ^iOi ^i^d find the recursion relation for the energy is: 

where label excited states of the unperturbed Hamitlonian. 

We can explicitly write out the first few terms of these series in terms of the unperturbed wave functions and 
energies. We drop all superscripts, which are 0. 

= (v^oli^ilV^o) 



E 



(2) _ 
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E,- 


Eo 


— ^ 


1 




i 




Eo 




1 





■(^o|Hi|V'o>(^,|i?i|V'o) 



E , - E, 



e^^^.j:(-^jI^\(^o\hm.)\ 



^ (^o|i^i|V'j)(^j|i?i|^.)(^.|^i|V'o) ) (74) 



{Ei — Eo){Ej — Eo 

The matrix elements are all ±t or 0, except for certain hoppings in the vicinity of the soliton, which have 

the form t{l + cos (/eg)). 

We are interested in corrections higher than 3'"'* order to investigate possible super-solid states, which consist of a 
background of filling and double occupancies at filling p/{nq). The non- vanishing energy corrections occur at 
multiples of n, as other orders cannot map the ground state back to itself. 

Here we consider the simplest case, n = 2. Since the first order correction to E vanishes, we may write i/''-^-* in the 
form 

l^^'^)-E rF ^ j,, \^^)(^^\Ht"\i'o) (75) 
j \Ej - Eq) 

where H^-^^ — e -Eq I^») ("0^1-^1 is the effective hopping Hamiltonian. Here \'4>i) has one occupied odd site. jV'j) 
differs from the ground state by 2 hoppings. There are two such terms: has occupation of even sites only, and 

\ipj-i) has two occupied odd sites. It is useful to separate these two terms explicitly: 

l^^'^)=E rF %> .>o)(^jo|gr^^lV'o)+^ ^ j,M n){^n\Hl^^\^o) (76) 
, {Ejo - Eq) ^ {Ej^ ~ Eq) 



The first piece looks like the first-order correction to a wave function on a re-scaled lattice, with an effective hopping 
coefficient given by H^^^ . The second piece does not contribute to the second-order energy correction. Hence in this 
notation, we also have 

=-(^,o|i/f//|V>o) (77) 
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(78) 



This is the leading-order energy correction. The third-order correction to the energy vanishes; the wave-function 
correction is: 

i ^ \ 2 

This gives the fourth order energy correction: 

j ^ 
= {i^o\H'ff\i'o){4'o\Hf^\^o)-J2 J, ^ ^ \{i'o\H'^"\i':.)\' 

\ mH'^"\i%)\' (79) 

where H'^^^ — J^i (^e -Eq)^ Hi\'4'i){''Pi\Hi. The last term is the analogue, in our effective theory, of a 'second-order' 
correction. The first two terms give an extra correction due to the structure of the underlying lattice. 



1. Matrix Elements and formulae used for numerics m the super-soUd states 



To evaluate the perturbative energy corrections, we must calculate the relevant matrix elements. Here we outline this 
calculation order by order for the super-solid phases, in which the results are somewhat complex. The corresponding 
calculations for the convex regime are a straightforward adaptation of the results of 0, and are not included here. To 
simplify the language, for the remainder of this section we will use qSS and CGS to refer to states with and without 
solitons in the doubly occupied sites, respectively. 

Second Order correction 

The second-order energy is given by 

E(^^ =-{i'o\Ht"\^Po) (80) 

The matrix elements 



(81) 



are \/2t when involves hopping from a doubly occupied site (as there are two identical particles to choose from), 
and t otherwise. (This is because 6i|no) = y/no\nQ — 1). A factor of V2 is also introduced by hopping onto an occupied 

site, since bl\nQ — 1) = y^?io|^o)- Hence hopping a particle off of a doubly occupied site, and then back on, incurs a 
factor of 2). 

Hence, for the CGS, 



E 



1 



Ei — Eq 



E 



Ei — Eq 



(82) 



where Ei corresponds to hopping the i*'* particle to the right. (The overall factor of 2 then accounts for the left 
hoppings). This can be written conveniently in the form: 



E 



(2) 



''i,R 



(83) 
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where t^'^ = e -Eq ^^^^ z — 1 is singly occupied, and e-^e„ doubly occupied. This accounts for the 

denominator when hopping a particle from site i — 1 one site to the right, to site i. We may equally well define f^// , 

which is ^.1^^ if site i + 1 is singly occupied, and ^.^^^ if it is doubly occupied. 

In the solitonized ground states, there is an additional contribution, in which the particle at one end of the soliton 
hops two sites to the right (or left) effectively hopping the soliton by 2q sites. This extra contribution is equivalent 
to the first order energy gain of the solitonized state in the non-super-solid phase. The t'^'^^ for this is the same as for 
this particle hopping out and returning to the same site. Thus for the qSS 

E% = -2 E ^" - 2^1^^ ^0^(2^9) (84) 

i 

where t^'^'^ contains the appropriate energy denominator for the hopping from site at the end of the soliton. The 
factor of 2 in front of the sum accounts for left and right hoppings; in the second term, the particle can hop from the 
right end of the soliton ending at site is to the left end of a (new) soliton beginning at ig -\- 2, or vice versa ( reversing 
left and right). The difference between the CGS and qSS energy corrections is thus: 

SE(^^ = -2 f t^^] - 2tt!f cos(2fcg) - ± (85) 

1=1 ^ i=l 

Here N is the number of particles on the chain. In practice the sum may be truncated at finite N as the difference 
between CGS and qSS contributions falls off rapidly far from a soliton. In this term, however, we have assumed that 
there is one particle in the CGS for each particle in the qSS; we must correct for this by subtracting the (average) 
energy of one particle in the CGS, given in the last term. 

Fourth order correction 

The fourth order energy correction is: 



E 



«0 



— ^|(Vo|ff^^^|^.o)P (86) 



The first term in the sum is a product of E^^\ calculated above, and E^-^^ given by 



{E, - EoY- ^^{E,- Eo) 

i^tDO iDO 



(87) 



Here again, the solitonized ground state has an extra contribution compared to the CGS, of the form ((M)) . Each term 
in this sum is uniquely labeled by choosing two site indices, i and j, to be the target sites of intermediate hopping, 
and a direction of hopping onto and off of each. 

In the second term, lipi) are states which arc connected to l^po) by hopping two particles k and j onto odd sites. 
The choice of k and j uniquely fixes the intermediate state IV'ii)- Note, however, that for each choice of k and j there 
are two possible hopping sequences: |i/'o) iV'fc) iV'i) IV'i) IV'o), and \tpo) \tpk) \ipk) \ipo), 

where \ipk) denotes the state in which particle k has hopped, but particle j has not. These two possible paths differ 
only in the order in which the hopping back to the ground state is performed. 

Hence the second term of ((86)) can be expressed as: 



1 1 



[Ei+i.j+i — EQ){Ej+i — i?o)2 (iJi+ij+i — Eo){Ei^i — Eo){Ej+i — Eq) 



E 



[Ei+i — Eo){Ej-^-i — EqY " [Ei-i^i+i — Eo){Ei+i — Eo){Ei^i — Eq) 



Here i and j are sites which are occupied by and Uj bosons in the ground state. -Ej+i are the energies of the 

intermediate states produced by hopping a boson off of sites i and j (by hopping to the right; the term for hopping to 
the left appears with a — sign). Ei^i,j^i represents the energy of the intermediate state with particles hopped from 
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both sites i and j. Here again, in the case of sohtonizcd states the particle at the end of the sohton can hop either 
'out and back', or 'out and out' ~ the latter producing a translation of the soliton. 

For both CGS and qSS states, the sum of the first two terms in Eq. ((86|) is small. In particular, the cross-terms 
between hoppings of the soliton itself and hoppings of the background lattice approximately cancel. This is as it 
should be - if they did not, these terms would give an energy splitting between qSS and CGS states which diverges 
in the thermodynamic limit! 

The final contribution resembles a second-order contribution with a re-scaled lattice constant and hopping terms 
which we will call e'^^ and E^^^ . The first of these has the form: 



e// 12 



Ei — Eq 



where i labels the locus of the hopped particle in the excited state, and i — 2 its locus in the ground state (prior to 
hopping). The factor of 2 accounts for hoppings to the left and right. The tl[_2 r is the effective hopping from site 
« — 2 to site i — 1 (which gives the same energy denominator as a leftward hop from site i to site i — 1). As in the 
non-super-solid case, there is one special intermediate configuration \ip*) in the solitonized case which permits two 
distinct hoppings which return to the ground state, as the ground state in this case is a momentum eigenstate. The 
intermediary configurations for the 'forward' and 'backward' hoppings arc mirror images of each other in this case. 
Hence for the solitonized ground state: 



ueff |2 ue//|2 



Ei~ Eq Ei- — Eq 

l^ls 



(89) 



We exclude hopping to the right from is in the first sum, as it merely hops the soliton, and has been included at 
second order. The special hopping from site i* gets a matrix element of 4, as it also occurs once in the first sum. 

The final contribution represents the effect of hopping both particles off of a doubly occupied site and then back, 
and has the form 



^ EiR,R — Eq Ei,RL — E{ 



{\t'i%? + nMi\) (90) 



where Ei is the energy of the intermediate configuration with both particles removed from the original site. Here one 
factor of 2 is, as usual, to account for both left and right hoppings. In the first term, both particles are hopped to 
the same site; the factor of 4 = (■\/2)'' arises because all hops move a particle either onto or off of a doubly occupied 
site. In the second term, one particle is hopped to the right and one to the left; consequently there is only a single 
factor of 2 due to double occupancies. However, the return hop may be executed in the same order as the outward 
hop (l^i'^P term), or the opposite order {{tl^^il^i \) ■ The form of this correction is the same for the CGS and for 
solitonized states. 



B. Bounds on the volume of Devil's staircase lost at small t 



Here we estimate the total volume of the liquid states in the phase diagram for small but finite t. 
First, let us review the situation at i = 0, where the total volume of liquid states is 0. The range of stability of a 
Mott lobe with denominator q is given by 

^nq{V{nq + l) + V{nq-l)-2V{nq)) (91) 

n 

To calculate the volume occupied by all rationally filled states, we sum over q. The multiplicity of each q is the 
number of rationals between and 1 with denominator q, which is given by the Totient function (l){q). This gives the 
volume: 



51 '^(9) ^h'^'^i'l'^ + 1) + qnV{qn - 1) - 2qnV{qn)] 

q n 



(92) 
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Switching the order of summation in the first term, and letting m = qn + I, gives: 



E 



g I rn — 1 



(m- l)V{m) 



(93) 



The other terms in (|92p can be treated similarly. The expression in square brackets is just (m — 1), so that the total 



IS 



^[(m - i f + {m+ i f - 2m^]V{m) = 2 ^ V{m) 



(94) 



which is in fact the threshold of stability (from below) of the integer-filled state. 

Now let us calculate the amount of this volume lost to superfluid states at small t. Since the form of the perturbative 
corrections is difficult to deduce exactly in q, we work to linear order in t. Hence from each Mott lobe of denominator 
q, a swath in chemical potential of length 2qt has been lost to superfluidity. Summing over q gives: 



Ai - 2t ^ q<f>iq) 
Q<Qo 



(95) 



To avoid subtracting off an infinite correction, we sum only over those values of q for which the Mott lobe has not, at 
this value of t, entirely disappeared. The threshold value of q is given by: 



2tqo = J2 



nq 



nq 



- 2- 



nq 



3g= 



{nq + 1)-^ [nq — If (nqf 
[q{2 + 3cot(7r/(7)2) - Stt cot(7r/g)(l + cot{Tr /qf)] 



(96) 



At t = 0; go is infinite; hence for small t go is very large, and we may expand the right-hand side to leading order in 
1/qo, giving: 



t 



15g5 



To complete the calculation, we must also estimate the volume of the lost Mott lobes: 



EE'^('?) 

q>qo n 

For go sufficiently large, we may use the bounds: 



nq 



nq 



nq 



(nq + If {nq — 1)'^ {nqf 



< (j){n) < n 



ln(n)^ 

and take the leading order in g, supposing go is large: 

nq ^ nq _ ^ nq 
{nq + If {nq — If {nqf 

Now we may use the bounds ([99]) to evaluate the sum over g: 

E'^(5)t774 < E 



E 



2tt^ 
15g4 



9>9o 



15g4 ^ 15g3 

q>qo 



The last sum can be evaluated exactly in terms of di-F functions; to leading order in 1/go, we obtain: 



A, < 



15go^ 



(97) 



(98) 



(99) 



(100) 



(101) 



(102) 



Calculating the lower bound is somewhat trickier; however, we may approximate it as an integral. The integral may 
be evaluated for integer values of the constant C, so we should round C up and then calculate. The result is: 



1 



2g2ln(go)C 



(103) 
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Note that here we do not expect the coefficient to be captured by the integral, only the qualitative behavior in go- 
Hence we have, for some constant a, 



< A2 < 



47r4 
15^ 



It remains to calculate Ai. To do so, we will first prove the identity: 



71 1 ^ 



fe=i 



k=l 



.Tl.o 3.71.9 



(104) 



(105) 



where fi{k) is the Mobius function, which takes on values of +1 for prime numbers, for perfect squares, and — 1 
otherwise. This can be shown by induction in n; the base case is true since 4>{1) ~ fJ.{l) = 1. For the inductive step, 
we must show: 

in + mn + 1) = i E Mfc) [l^j^ - ItJ' + ' L^J^ - 'ItJ' 



k=l 



1 I rt + 1 1 I « 



+ (n + l)/i(n + 1) 



^ow, [-^^^J" — Lf J° vanishes except when fc|n + 1, in which case it gives 



n+1 



n + 1 



n + 1 
k 



o" + l . .." + 1 3 1 



Substituting this into the series above, we obtain: 
(n + l)0(n + l) = i 

k\n+l,k<n+l 

+ (n + l)/i(n + 1) 

k\n 

= (n + l)0(n + l) 

where the last equality is a basic identity of Totient functions. (Proved on Wikipedia). 
Now we may use the fact that 



n 
k 



, n , n 
'<^-k^^-k 



to put bounds on the series. The upper bound is: 



tEt^ + YE^+fiE^(^) 
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while the lower bound is: 



n"^ \ ^ u(k) v? uik) n , 
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To estimate these contributions, we use the relations: 



f 4^1 8+0(l/„) 
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Hence at large n, we have 



km = i-^+O (login) /n)) (113) 



k=l 

Plugging this into the expression for Ai above, we have 
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(114) 



Plugging in for t in terms of qo, and keeping only the leading term, we obtain: 



2 



47r 

Adding up both contributions, wc obtain the following asymptotic form for the volume of the liquid states in fi: 



ggln(qo)'^ 15g^ V 9o / ISgg 15g^ V 9o 

We can see that we expect the result to scale approximately as ^, or as t^/s. 

90 



C. Creation of double occupancies in CGS states: calculation of dipolar and quadrupolar interactions 

Here we discuss the formation of double occupancies in the CGS. We will make use of the fact that a double 
occupancy can always be formed by pinching about a hole in the CGS - that is, by moving some number of particles 
on the right of that hole leftwards, and some number of particles to the left of that hole rightwards. 

Recall that by even states, we mean states in which either p or g is even; states for which both are odd are odd 
states. In an even state, the pinching is symmetrical; in an odd state it cannot be. This produces a quadrupolar 
interaction between DO in even states, but a dipolar interaction for odd states. 

The pinching operation described in Sect. IVI A2l can be described by: 

Ri ^ Ri — cLi 

Lj Lj - bj (117) 

where i?^, Lj are the positions of the particles relative to the hole in the initial configuration from which the double 
occupancy is formed. (For an odd state, Ri has one more element than Li.) 

Now let us consider interactions between double occupancies. The interaction energy between two DO formed a 
distance d apart is given by the difference in energies of the configuration before and after all distances Ri , Li have 
been altered, minus the difference in energy when only a single DO has been formed. The relevant terms are the 
interaction energies between particles that have been pinched to form the first defect with particles that have been 
pinched to form the second: 

^ - [l/{d - R, - + 6,)' + ^/{d - R. - L, + a,)3] 

- l/{d - R, - Lj +a, + hjf -l/{d- R, - L^)^] 

- [l/{d + R, + Lj + b,)^ + l/{d + R^ + Lj + aj)^ 

- l/{d + R, + Lj + aj + b,)^ - l/{d + R, + Lj)^ 

+ [l/{d + i?j - Rj + bj - 6,)3 - l/{d + i?j - Rj - 6,)3 

- l/{d + R,- Rj + bjf + l/{d + R, - Rj)^] 

+ [l/{d + L, - Lj + aj - a,)'^ - l/{d + L, - Lj - a,)^ 

- l/{d + L, - Lj + aj)^ + l/{d + L, - L,)3] (118) 
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For the even case, Ri = Li and ai = bi, giving 

J2 - [^/{d - Rij + aif + - R^J + ajf - l/{d - i?y + a, + a,f - l/{d - R,,f] 

- [l/{d + i?y + a^f + l/{d + + a^f - l/{d + R,jf - l/{d + R,^ + a, + aj^] 
+ 2[l/(d + - 2Lj + aj - a,)3 - + R,^ - 2Lj - aj)^ 

- + i?y - 2Lj + ajf + l/{d + R,j - 2Ljf] (119) 

where Rij = Li + Lj. Each term in square brackets is negative definite, by convexity. A straightforward manipulation 
of the standard definition of convexity shows that for a > b and any x, we have 

V{x + a) + V{x + b)> V{x + c) + V{x + {a + b-c)) . (120) 

Equally, for interactions, the last square bracketed term is necessarily less than the sum of the first two (since 
these lie 'closer in', if you will) and hence the overall interaction is repulsive. 

For the odd case, again all terms in square brackets are negative (by convexity). However, as there are more Ri 
than Li, there are more terms in the last 2 lines than in the first two. (One more term, to be precise). Since the last 
2 lines are negative, this results in an interaction that is repulsive at long distances. 

At short distances, however, the interaction is dominated by the first line, which gives a positive contribution (and 
hence the double occupancies repel at sufficiently short distance scales). Certainly, when max{Ri + Lj) « d, we expect 
the potential to be repulsive and the minimum therefore occurs at some d > da = max{Ri + Lj). The transition is 
thus to a density of double occupancies that sit at least do sites apart, and the transition is to a density of double 
occupancies that is less than l/do- 

Further, do should increase linearly in q. This is easiest to see by means of example, but basically creating a DO 
will form hole-like solitons, by which I mean you can change distances of q to g + 1, but not to q + 2. To arrive at 
such a configuration involves a re-arrangement over q sites to the right of the DO (and something like q/2 to the left, 
seemingly) . 

As an example, we plot the 2 DO interaction potential as a function of separation between the doubly occupied 
site of each DO for the 1/3 and 1/5 filled states in Fig. [51 The minima sit at 15 and 75 lattice spacings, respectively. 

Hence as Uq is lowered, we expect even states to undergo a second-order phase transition, as the interaction between 
DO is quadrupolar and hence repulsive. For odd states the interaction is attractive, suggesting a first-order phase 
transition in which a finite density of defects forms. However, this density of defects vanishes at least as l/q, and 
quite possibly more rapidly than this. 



D. Lattice-scale arguments for adding charge as double occupancies 

Here we present a more complete argument that charges added as solitons have a stronger mutual repulsion en- 
ergy than charges added as double occupancies. At long length scales we may employ a simple intuition based on 
electrostatics. That is, when adding charges to the system we must compare both their self-energy (in this case, the 
energy of adding a single charge) and their interaction. The interaction energy between n charges q, spread over some 
finite interval of length L in the system, is E{n) q (^) . This is always less than the interaction energy between mn 

charges of charge q/m for any integer to > 1, which is E{m, n) ^ ( ) ^ = m^E{n). We work in the regime where 
the self-energy favors doubly occupied states over solitons - hence the charge will never fractionalize. 

Here we outline a lattice argument to show that this intuition gives the correct result even at short distances. We 
treat only the case i' > 1/2, as the form of the doubly occupied states is simpler to describe; however, we expect 
similar arguments to hold for arbitrary filling fractions. 

Consider adding a finite density of charges to a state with existing double-occupancies. We break the energy cost 
of introducing the new particle into two pieces: 

E = Ei + Ei (121) 

where Ei is the self-energy, and Ei is the interaction energy. The interaction energy contains the repulsion between the 
extra charge and any existing double-occupancies. If the extra charge breaks into solitons, Ej also contains repulsion 
terms between these solitons. Since we add a finite density of charge, we can consider that each charge is effectively 
confined to a finite region of length r in the system. 

Consider first the 1/2-filled state. If the charge fractionalizcs, two solitons will be formed. We place these solitons at 
ri and r, respectively, from the double occupancy. Here ri and r measure the distance from the DO to the outer-most 
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FIG. 8. Two-DO interaction potential Vd for the 1/3 (lower plot) and 1/5 (upper plot) filled states. Here roo is the distance 
between the two double occupancies in units of the lattice constant. The minima of these potentials are at roo ~ 75 and 15, 
respectively, for the 1/3 and 1/5 -filled states. 



edge of each soliton. One soliton sits as far from the DO as possible, as the repulsion between the DO and each soliton 
is greater than the soliton-soliton repulsion, as a DO has integral charge while the charge of a soliton is fractionalized. 
We may think of forming this arrangement by excising a hole at ri, pushing the rest of the pattern inwards, and 
placing an extra particle at r. For r = 10, this looks like: 



20101010101 
20101101010 
20101101011 

Note that ri is necessarily odd, and r here must be even. Then clearly we have: 

(r-ri-l)/2 



(122) 
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(123) 

(r - n + 2ri + 1)3 (r-n + 2n + 2)3j ^ ^ 

The first line here is the interaction between the soliton and the DO; the second line gives the soliton-soliton repulsion. 
(One pair of terms for each particle in the soliton) . The sums are all positive, since the l/r^ potential is monotonically 
decreasing. Hence, the interaction energy from adding a particle as solitons is strictly greater than that for adding it 
as a double occupancy: 



Ei{sol) > 



1 



El {DO) 



(124) 



Since r above is arbitrary, this argument also applies to the case of adding a particle to a system with multiple 
double occupancies in the 1/2-filled state. Hence in any finite system, we can induct on the number of particles to 
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show that every fihing p/q> 1/2 consists only of double occupancies in the 1/2 filled state, and never of soliton-like 
insertions ...0110.... In an infinite system, an infinite number of particles must be added at once in order to change the 
filling fraction; in this case the particles want to spread out homogeneously and hence we again confine them to within 
some distance r of existing double occupancies. (Similar arguments show that the repulsion between four solitons in 
a finite region is greater than the repulsion between two double-occupancies, as the latter can spread farther apart). 

Finally, we may generalize this argument to arbitrary fillings p/q > 1/2. Again we form solitons by pushing a 
hole from some radius ri to the farthest possible radius r from the DO, and inserting an extra particle at r. This is 
illustrated below for the 2/3 filled state: 

2011011011011011 
2011011110110110 

2011011101110111 (125) 

We may hop some (but not all) of the particles back to their original positions to form solitons, as shown in the third 
line. The energy of such a state clearly obeys 

Ei{sol) > ^ = Ei{DO) (126) 

and hence again, all additional particles will enter as double occupancies. This shows that the doubly occupied state 
is locally stable at all length scales. 
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